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In this review, we describe and analyze a mesoscale simulation method for fluid 
flow, which was introduced by Malevanets and Kapral in 1999, and is now called 
multi-particlo collision dynamics (MFC) or stochastic rotation dynamics (SRD). The 
method consists of alternating streaming and collision steps in an ensemble of point 
particles. The multi-particle collisions are performed by grouping particles in colli- 
sion cells, and mass, momentum, and energy are locally conserved. This simulation 
technique captures both full hydrodynamic interactions and thermal fluctuations. 
The first part of the review begins with a description of several widely used MFC 
algorithms and then discusses important features of the original SRD algorithm and 
frequently used variations. Two complementary approaches for deriving the hydro- 
dynamic equations and evaluating the transport coefficients are reviewed. It is then 
shown how MFC algorithms can be generalized to model non-ideal fluids, and binary 
mixtures with a consolute point. The importance of angular-momentum conserva- 
tion for systems like phase-separated liquids with different viscosities is discussed. 
The second part of the review describes a number of recent applications of MFC 
algorithms to study colloid and polymer dynamics, the behavior of vesicles and cells 
in hydrodynamic flows, and the dynamics of viscoelastic fluids. 

PACS number(s): 47.ll.-j, 05.40.-a, 02.70.Ns 

1 

Introduction 

"Soft Matter" is a relatively new field of research that encompasses tradi- 
tional complex fluids such as amphiphilic mixtures, colloidal suspensions, and 
polymer solutions, as well as a wide range of phenomena including chemically 
reactive flows (combustion), the fluid dynamics of self-propelled objects, and 
the visco-elastic behavior of networks in cells. One characteristic feature of 
all these systems is that phenomena of interest typically occur on mesoscopic 
length-scales ranging from nano- to micrometers — and at energy scales com- 
parable to the thermal energy ksT. 
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Because of the complexity of these systems, simulations have played a par- 
ticularly important role in soft matter research. These systems are challenging 
for conventional simulation techniques due to the presence of disparate time, 
length, and energy scales. Biological systems present additional challenges be- 
cause they are often far from equilibrium and are driven by strong spatially 
and temporally varying forces. The modeling of these systems often requires 
the use of "coarse-grained" or mesoscopic approaches that mimic the behavior 
of atomistic systems on the length scales of interest. The goal is to incorporate 
the essential features of the microscopic physics in models which arc compu- 
tationally efficient and are easily implemented in complex geometries and on 
parallel computers, and can be used to predict behavior, test physical the- 
ories, and provide feedback for the design and analysis of experiments and 
industrial applications. 

In many situations, a simple continuum description based on the Navier- 
Stokcs equation is not sufficient, since molecular- level details — including ther- 
mal fluctuations — play a central role in determining the dynamic behavior. 
A key issue is to resolve the interplay between thermal fluctuations, hydro- 
dynamic interactions, and spatio-temporally varying forces. One well-known 
example of such systems are microemulsions — a dynamic bicontinuous net- 
work of intertwined mesoscopic patches of oil and water — where thermal fluc- 
tuations play a central role in creating this phase. Other examples include 
flexible polymers in solution, where the coil state and stretching elasticity are 
due to the large configurational entropy. On the other hand, atomistic molec- 
ular dynamics simulations retain too many microscopic degrees of freedom, 
consequently requiring very small time steps in order to resolve the high fre- 
quency modes. This makes it impossible to study long timescale behavior such 
as self-assembly and other mcsoscalc phenomena. 

In order to overcome these difficulties, considerable effort has been devoted 
to the development of mesoscale simulation methods such as Dissipative Par- 
ticle Dynamics [1 3], Lattice-Boltzmann [4 6], and Direct Simulation Monte 
Carlo [7-9]. The common approach of all these methods is to "average out" 
irrelevant microscopic details in order to achieve high computational efficiency 
while keeping the essential features of the microscopic physics on the length 
scales of interest. Applying these ideas to suspensions leads to a simplifled, 
coarse-grained description of the solvent degrees of freedom, in which embed- 
ded macromolcculcs such as polymers axe treated by conventional molecular 
dynamics simulations. 

All these approaches are essentially alternative ways of solving the Navier- 
Stokes equation and its generalizations. This is because the hydrodynamic 
equations are expressions for the local conservation laws of mass, momentum, 
and energy, complemented by constitutive relations which reflect some aspects 
of the microscopic details. Frisch et al. [10] demonstrated that discrete algo- 
rithms can be constructed which recover the Navier-Stokes equation in the 
continuum limit as long as these conservation laws are obeyed and space is 
discretized in a sufficiently symmetric manner. 
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The first model of this type was a cellular automaton, called the Lattice- 
Gas- Automaton (LG). The algorithm consists of particles which jump between 
nodes of a regular lattice at discrete time intervals. Collisions occur when 
more than one particle jumps to the same node, and collision rules are cho- 
sen which impose mass and momentum conservation. The Lattice-Boltzmann 
method (LB) which follows the evolution of the single-particle probability 
distribution at each node — was a natural generalization of this approach. LB 
solves the Boltzmann equation on a lattice with a small set of discrete veloci- 
ties determined by the lattice structure. The price for obtaining this efficiency 
is numerical instabilities in certain parameter ranges. Furthermore, as origi- 
nally formulated, LB did not contain any thermal fluctuations. It became clear 
only very recently (and only for simple liquids) how to restore fluctuations by 
introducing additional noise terms to the algorithm [11]. 

Except for conservation laws and symmetry requirements, there are rel- 
atively few constraints on the structure of mesoscale algorithms. However, 
the constitutive relations and the transport coefficients depend on the details 
of the algorithm, so that the temperature and density dependencies of the 
transport coefficients can be quite different from those of real gases or liq- 
uids. However, this is not a problem as long as the functional form of the 
resulting hydrodynamic equations is correct. The mapping to real systems is 
achieved by tuning the relevant characteristic numbers, such as the Reynolds 
and Peclet numbers [12, 13], to those of a given experiment. When it is not 
possible to match all characteristic numbers, one concentrates on those which 
are of order unity, since this indicates that there is a delicate balance between 
two effects which need to be reproduced by the simulation. On occasion, this 
can be difficult, since changing one internal parameter, such as the mean free 
path, usually affects all transport coefficients in different ways, and it may 
happen that a given mesoscale algorithm is not at all suited for a given ap- 
plication [14-17]. 

In this review we focus on the development and application of a particle- 
based mesoscopic simulation technique which was recently introduced by Mal- 
evanets and Kapral [18,19]. The algorithm, which consists of discrete stream- 
ing and collision steps, shares many features with Bird's Direct Simulation 
Monte Garlo (DSMG) approach [7]. Gollisions occur at fixed discrete time in- 
tervals, and although space is discretized into cells to define the multi-particle 
collision environment, both the spatial coordinates and the velocities of the 
particles are continuous variables. Because of this, the algorithm exhibits un- 
conditional numerical stability and has an i?-theorem [18,20]. In this review, 
we will use the name multi-particle collision dynamics (MFC) to refer to this 
class of algorithms. In the original and most widely used version of MPG, 
collisions consist of a stochastic rotation of the relative velocities of the par- 
ticles in a collision cell. We will refer to this algorithm as stochastic rotation 
dynamics (SRD) in the following. 

One important feature of MPG algorithms is that the dynamics is well- 
defined for an arbitrary time step. At. In contrast to methods such as molecu- 
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lar dynamics simulations (MD) or dissipative particle dynamics (DPD) , which 
approximate the continuous-time dynamics of a system, the time step does 
not have to be small. MPC defines a discrctc-timc dynamics which has been 
shown to yield the correct long-time hydrodynamics; one consequence of the 
discrete dynamics is that the transport coefficients depend explicitly on At. In 
fact, this freedom can be used to tunc the Schmidt number, Sc [15]; keeping 
all other parameters fixed, decreasing At leads to an increase in Sc. For small 
time steps, Sc is larger than unity (as in a dense fluid), while for large time 
steps. Sc is of order unity, as in a gas. 

Because of its simplicity, SRD can be considered an "Ising model" for hy- 
drodynamics, since it is Galilean invariant (when a random grid shift of the 
collision cells is performed before each collision step [21]) and incorporates all 
the essential dynamical properties in an algorithm which is remarkably easy 
to analyze. In addition to the conservation of momentum and mass, SRD also 
locally conserves energy, which enables simulations in the microcanonical en- 
semble. It also fully incorporates both thermal fluctuations and hydrodynamic 
interactions. Other more established methods, such as Brownian Dynamics 
(BD) can also be augmented to include hydrodynamic interactions. However, 
the additional computational costs are often prohibitive [22,23]. In addition, 
hydrodynamic interactions can be easily switched off in MPC algorithms, 
making it easy to study the importance of hydrodynamic interactions [24,25]. 

It must, however, be emphasized that all local algorithms such as MPC, 
DPD, and LB model compressible fluids, so that it takes time for the hydro- 
dynamic interactions to "propagate" over longer distances. As a consequence, 
these methods become quite inefficient in the Stokes limit, where the Reynolds 
number approaches zero. Algorithms which incorporate an Oseen tensor do 
not share this shortcoming. 

The simplicity of the SRD algorithm has made it possible to derive ana- 
lytic expressions for the transport coefficients which are valid for both large 
and small mean free paths [26-28] . This is usually very difficult to do for other 
mesoscale particle-based algorithms. Take DPD as an example: the viscosity 
measured in Ref. [29] is about 50% smaller than the value predicted theo- 
retically in the same paper. For SRD, the agreement is generally better than 
1%. 

MPC is particularly well suited for (i) studying phenomena where both 

thermal fluctuations and hydrodynamics arc important, (ii) for systems with 
Reynolds and Peclet numbers of order 0.1 to 10, (iii) if exact analytical expres- 
sions for the transport coefficients and consistent thermodynamics are needed, 
and (iv) for modeling complex phenomena for which the constitutive relations 
are not known. Examples include chemically reacting flows, self-propelled ob- 
jects, or solutions with embedded macromolecules and aggregates. 

If thermal fluctuations arc not essential or undesirable, a more traditional 
method such as a flnite-element solver or a Lattice-Boltzmann approach is 
recommended. If, on the other hand, inertia and fully resolved hydrodynamics 
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are not crucial, but fluctuations are, one might be better served using Langevin 
or Brownian Dynamics. 

This review consists of two parts. The first part begins with a description 
of several widely used MPC algorithms in Sec. 2, and then discusses impor- 
tant features of the original SRD algorithm and a frequently used variation 
(MPC- AT) , which effectively thermostats the system by replacing the relative 
velocities of particles in a collision cell with newly generated Gaussian random 
numbers in the collision step. After a qualitative discussion of the static and 
dynamic properties of MPC fluids in Sec. 3, two alternative approaches for 
deriving the hydrodynamic equations and evaluating the transport coefficients 
are described. First, in Sec. 4, discrete-time projection operator methods are 
discussed and the explicit form of the resulting Creen-Kubo relations for the 
transport coefficients are given and evaluated. Subsequently, in Sec. 5, an 
alternative non-equilibrium approach is described. The two approaches com- 
plement each other, and the predictions of both methods arc shown to be in 
complete agreement. It is then shown in Sec. 6 how MPC algorithms can be 
generalized to model non-ideal fluids and binary mixtures. Finally, various 
approaches for implementing slip and no-slip boundary conditions — as well as 
the coupling of embedded objects to a MPC solvent — are described in Sec. 7. 
In Sec. 8, the importance of angular-momentum conservation is discussed, 
in particular in systems of phase-separated fluids with different viscosities 
under flow. An important aspect of mesoscale simulations is the possibility 
to directly assert the effect of hydrodynamic interactions by switching them 
off, while retaining the same thermal fluctuations and similar friction coeffi- 
cients; in MPC, this can be done very efficiently by an algorithm described 
in Sec. 9. The second part of the review describes a number of recent appli- 
cations of MPC algorithms to study colloid and polymer dynamics, and the 
behavior of vesicles and cells in hydrodynamic flows. Sec. 10 focuses on the 
non-equilibrium behavior of colloidal suspensions, the dynamics of dilute so- 
lutions of linear polymers both in equilibrium and under flow conditions, and 
the properties of star polymers — also called ultra-soft colloids — in shear flow. 
Sec. 11 is devoted to the review of recent simulation results for membranes in 
flow. After a short introduction to the modeling of membranes with different 
levels of coarse-graining, the behavior of fluid vesicles and red blood cells, 
both in shear and capillary flow, is discussed. Finally, a simple extension of 
MPC for viscoclastic solvents is described in Sec. 12, where the point particles 
of MPC for Newtonian fluids axe replaced by harmonic dumbbells. 

A discussion of several complementary applications — such as chemically 
reactive flows and self-propelled objects — can be found in a recent review of 
MPC by R. Kapral [30]. 
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2 

Algorithms 

In the following, we use the term multi-particle collision dynamics (MPC) to 
describe the generic class of particle-based algorithms for fluid flow which con- 
sist of successive free-streaming and multi-particle collision steps. The name 
stochastic rotation dynamics (SRD) is reserved for the most widely used algo- 
rithm which was introduced by Malevanets and Kapral [18]. The name refers 
to the fact that the collisions consist of a random rotaf/ion of the relative ve- 
locities 6vi = Vj — u of the particles in a collision cell, where u is the mean 
velocity of all particles in a cell. There are a number of other MPC algorithms 
with different collision rules [31-33]. For example, one class of algorithms uses 
modified collision rules which provide a nontrivial "collisional" contribution 
to the equation of state [33,34]. As a result, these models can be used to model 
non-ideal fluids or multi-component mixtures with a consolute point. 



Stochastic Rotation Dynamics (SRD) 

In SRD, the solvent is modeled by a large number A'' of point-hke particles 
of mass m which move in continuous space with a continuous distribution of 
velocities. The algorithm consists of individual streaming and collision steps. 
In the streaming step, the coordinates, ri{t), of all solvent particles at time t 
are simultaneously updated according to 



where Vi(t) is the velocity of particle i at time t and At is the value of the 
discretized time step. 

In order to define the colUsions, particles are sorted into cells, and they 
interact only with members of their own cell. Typically, the system is coarse- 
grained into cells of a regular, typically cubic, grid with lattice constant a. 
In practice, lengths are often measured in units of o, which corresponds to 
setting a = 1. The average number of particles per cell, M, is typically chosen 
to be between three and 20. The actual number of particles in cell at a given 
time, which fluctuates, will be denoted by Nc. The collision step consists of a 
random rotation R of the relative velocities ^Vj = Vj — u of all the particles 
in the collision cell. 



All particles in the cell are subject to the same rotation, but the rotations 
in different cells and at different times are statistically independent. There 
is a great deal of frcxxloni in how the rotation step is implemented, and any 
stochastic rotation matrix which satisfies semi-detailed balance can be used. 
Here, we describe the most commonly used algorithm. In two dimensions, R 
is a rotation by an angle ±a, with probability 1/2. In three dimensions, a 



2.1 



ri{t + At)=ri{t)+AtVi{t), 



(1) 



Vi{t + At) = u{t) + R • 5vi{t) . 



(2) 
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rotation by a fixed angle a about a randomly chosen axis is typically used. 
Note that rotations by an angle —a need not be considered, since this amounts 
to a rotation by an angle a about an axis with the opposite orientation. If we 
denote the randomly chosen rotation axis by R, the explicit collision rule in 
three dimensions is 

v,{t + At) u{t) + Svi^j_{t) cos(a) 

+ {Svi,± {t) X R) sin(a) + <5vi, y (i) , (3) 

where _L and || are the components of the vector which are perpendicular 
and parallel to the random axis R, respectively. Malevanets and Kapral [18] 
have shown that there is an ff-theorem for the algorithm, that the equilib- 
rium distribution of velocities is Maxwellian, and that it yields the correct 
hydrodynamic equations with an ideal-gas equation of state. 

In its original form [18, 19], the SRD algorithm was not Galilean invariant. 
This is most pronounced at low temperatures or small time steps, where the 
mean free path, A = At^JksT /m. is smaller than the cell size a. If the particles 
travel a distance between collisions which is small compared to the cell size, 
essentially the same particles collide repeatedly before other particles enter 
the cell or some of the participating particles leave the cell. For small A, large 
numbers of particles in a given cell remain correlated over several time steps. 
This leads to a breakdown of the molecular chaos assumption — i.e., particles 
become correlated and retain information of previous encounters. Since these 
correlations are changed by a homogeneous imposed flow field, V, Galilean 
invariance is destroyed, and the transport coefficients depend on both the 
magnitude and direction of V. 

Ihle and KroU [20, 21] showed that Galilean invariance can be restored 
by performing a random shift of the entire computational grid before every 
collision step. The grid shift constantly groups particles into new collision 
neighborhoods; the collision environment no longer depends on the magni- 
tude of an imposed homogeneous flow fleld, and the resulting hydrodynamic 
equations are Galilean invariant for arbitrary temperatures and A/[ach number. 
This procedure is implemented by shifting all particles by the same random 
vector with components uniformly distributed in the interval [—a/2, a/2] be- 
fore the collision step. Particles are then shifted back to their original positions 
after the collision. 

In addition to restoring Galilean invariance, this grid-shift procedure accel- 
erates momentum transfer between cells and leads to a coUisional contribution 
to the transport coefficients. If the mean free path A is larger than a/2, the 
violation of Galilean invariance without grid shift is negligible, and it is not 
necessary to use this procedure. 
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2.1.1 

SRD with Angular Momentum Conservation 

As noted by Pooley and Yeomans [35] and confirmed in Ref. [28] , the macro- 
scopic stress tensor of SRD is not symmetric in daV/s- The reason for this 
is that the multi-particle collisions do not, in general, conserve angular mo- 
mentum. The problem is particularly pronounced for small mean free paths, 
where asymmetric coUisional contributions to the stress tensor doniinate the 
viscosity (see Sec. 4.1.1). In contrast, for mean free paths larger than the cell 
size, where kinetic contributions dominate, the effect is negligible. 

An anisotropic stress tensor means that there is non-zero dissipation if 
the entire fluid undergoes a rigid-body rotation, which is clearly unphysical. 
However, as emphasized in Ref. [28], this asymmetry is not a problem for 
most applications in the incompressible (or small Mach number) limit, since 
the form of the Navier-Stokes equation is not changed. This is in accordance 
with results obtained in SRD simulations of vortex shedding behind an obsta- 
cle [36] , and vesicle [37] and polymer dynamics [14] . In particular, it has been 
shown that the linearized hydrodynamic modes are completely unaffected in 
two dimensions; in three dimensions only the sound damping is slightly mod- 
ified [28]. 

However, very recently Gotze et al. [38] identified several situations involv- 
ing rotating flow fields in which this asymmetry leads to significant deviations 
from the behavior of a Newtonian fluid. This includes (i) systems in which 
boundary conditions are defined by torques rather than prescribed velocities, 
(ii) mixtures of liquids with a viscosity contrast, and (iii) polymers with a 
locally high monomer density and a monomer-monomer distance on the order 
of or smaller than the lattice constant, a, embedded in a MFC fluid. A more 
detailed discussion will be presented in Sec. 8 below. 

For the SRD algorithm, it is possible to restore angular momentum con- 
servation by having the collision angle depend on the specific positions of the 
particles within a collision cell. Such a modification was first suggested by Ry- 
der [39] for SRD in two dimensions. She showed that the angular momentum 
of the particles in a collision cell is conserved if the collision angle a is chosen 
such that 

sin(a) = -2AB/(A2 + B^) and cos(a) = {A^ - B^) /{A^ + B^) (4) 
where 

A = ^[rj X (vi -u)]|^ and B = ^ • (v^ - u). (5) 
1 1 

When the collision angles are determined in this way, the viscous stress tensor 
is symmetric. Note, however, that evaluating Eq. (4) is time-consuming, since 
the collision angle needs to be computed for every collision cell every time 
step. This typically increases the CFU time by a factor close to two. 
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A general procedure for implementing angular-momentum conservation in 
multi-particle collision algorithms was introduced by Noguchi et al. [32]; it is 
discussed in the following section. 

2.2 

Multi-Particle Collision Dynamics with Anderson Thermostat (MPC-AT) 

A stochastic rotation of the particle velocities relative to the center-of-mass 
velocity is not the only possibility for performing multi-particle collisions. 
In particular, MPC simulations can be performed directly in the canonical 
ensemble by employing an Anderson thermostat (AT) [31,32]; the resulting 
algorithm will be referred to as MPC-AT-a. In this algorithm, instead of per- 
forming a rotation of the relative velocities, {^v,}, in the collision step, new 
relative velocities are generated. The components of {i5v™"} are Gaussian 
random numbers with variance -^/fesT/m. The collision rule is [32,38] 

Vi{t + At) = u{t) + 5v,™" = u(i) + v^" - J2 v™"/iVc , (6) 

jGcell 

where A^c is the number of particles in the collision cell, and the sum runs over 
all particles in the cell. It is important to note that MPC-AT is both a colli- 
sion procedure and a thermostat. Simulations are performed in the canonical 
ensemble, and no additional velocity rescaling is required in non-equilibrium 
simulations, where there is viscous heating. 

Just as SRD, this algorithm conserves momentum at the cell level but 
not angular momentum. Angular momentum conservation can be restored 
[32, 39] by imposing constraints on the new relative velocities. This leads to 
an angular-momentum conserving modification of MPC-AT [32,38], denoted 
MPC-AT+a. The collision rule in this case is 

Vi{t + At) = U{t) + Vi^ran " ^ '^i,ranlNc 

cell 

+ \ mn-' ^ [r,,e X (v, - V™")] X r,,, i , (7) 

{ jecell J 

where JT is the moment of inertia tensor of the particles in the cell, and 

Ti.c = ^ R-c is the relative position of particle i in the cell and Rc is the 
center of mass of all particles in the cell. 

When implementing this algorithm, an unbiased multi-particle collision is 
first performed, which typically leads to a small change of angular momen- 
tum, AL. By solving the linear equation — AL = Tl u3, the angular velocity u> 
which is needed to cancel the initial change of angular momentum is then de- 
termined. The last term in Eq. (7) restores this angular momentum deficiency. 
MPC-AT can be adapted for simulations in the micro-canonical ensemble by 
imposing an additional constraint on the values of the new random relative 
velocities [32]. 
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2.2.1 

Comparison of SRD and MPC-AT 

Because d Gaussian random numbers per particle are required at every it- 
eration, where d is the spatial dimension, the speed of the random number 
generator is the limiting factor for MPC-AT. In contrast, the efficiency of 
SRD is rather insensitive to the speed of the random number generator since 
only d — I uniformly distributed random numbers arc needed in every box 
per iteration, and even a low quality random number generator is sufficient, 
because the dynamics is self-averaging. A comparison for two-dimensional sys- 
tems shows that MPC-AT-a is about a factor 2 to 3 times slower than SRD, 
and that MPC-AT+a is about a factor 1.3 to 1.5 slower than MPC-AT-a [40]. 

One important difference between SRD and MPC-AT is the fact that re- 
laxation times in MPC-AT generally decrease when the number of particles 
per cell is increased, while they increase for SRD. A longer relaxation time 
means that a larger number of time steps is required for transport coefficients 
to reach their asymptotic value. This could be of importance when fast oscil- 
latory or transient processes are investigated. As a consequence, when using 
SRD, the average number of particles per cell should be in range of 3 — 20; 
otherwise, the internal relaxation times could be no longer negligible com- 
pared to physical time scales. No such limitation exists for MPC-AT, where 
the relaxation times scale as (lnM)~^, where M is the average number of 
particles in a collision cell. 

2.3 

Computationally Efficient Cell-Level Thermostating for SRD 

The MPC-AT algorithm discussed in Sec. 2.2 provides a very efEc;ient particle- 
level thermostating of the system. However, it is considerably slower than 
the original SRD algorithm, and there are situations in which the additional 
freedom offered by the choice of SRD collision angle can be useful. 

Thermostating is required in any non-equilibrium MPC simulation, where 
there is viscous heating. A basic requirement of any thermostat is that it 
does not violate local momentum conservation, smear out local flow profiles, 
or distort the velocity distribution too much. When there is homogeneous 
heating, the simplest way to maintain a constant temperature is to just rescale 
velocity components by a scale factor S, u"*^™ = Sv^, which adjusts the total 
kinetic energy to the desired value. This can be done with just a single global 
scale factor, or a local factor which is different in every cell. For a known 
macroscopic flow profile, u, like in shear flow, the relative velocities v — u can 
be rescaled. This is known as a profile-unbiased thermostat; however, it has 
been shown to have deficiencies in molecular dynamics simulations [41]. 

Here we describe an alternative thermostat which exactly conserves mo- 
mentum in every cell and is easily incorporated into the MPC collision step. It 
was originally developed by Heyes for constant-temperature molecular dynam- 
ics simulations; however, the original algorithm described in Ref. [42] violates 
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detailed balance. The thermostat consists of the following procedure which is 
performed independently in every collision cell as part of the collision step. 

1. Randomly select a real number tp G [1, 1 + c], where c is a small number 
between 0.05 and 0.3 which determines the strength of the thermostat. 

2. Accept this number as a scaling factor S = ip with probability 1/2; oth- 
erwise, take S = 1/tp. 

3. Create another random number ^ S [0,1]. Rcscalc the velocities if ^ is 
smaller than the acceptance probability pa = min(l, A), where 



d is the spatial dimension, and Nc is the number of particles in the cell. 
The prcfactor in Eq. (8) is an cntropic contribution which accounts for 
the fact that the phase-space volume changes if the velocities are rescaled. 
4. If the attempt is accepted, perform a stochastic rotation with the scaled 
rotation matrix S R. Otherwise, use the rotation matrix R. 

This thermostat reproduces the Maxwell velocity distribution and docs not 
change the viscosity of the fluid. It gives excellent equilibration, and the devi- 
ation of the measured kinetic temperature from To is smaller than 0.01%. The 
parameter c controls the rate at whic;h the kinetic temperature relaxes to To, 
and in agreement with experience from MC-simulations, an acceptance rate in 
the range of 50% to 65% leads to the fastest relaxation. For these acceptance 
rates, the relaxation time is of order 5—10 time steps. The c;orresponding 
value for c depends on the particle number N^, in two dimensions, it is about 
0.3 for Nc = 7 and decreases to 0.05 for N^. = 100. This thermostat has been 
successfully applied to SRD simulations of sedimenting charged colloids [16]. 



Qualitative Discussion of Static and Dynamic Properties 

The previous section outlines several multi-particle algorithms. A detailed 
discussion of the link between the microscopic dynamics described by Eqs. (1) 
and (2) or (3) and the macroscopic hydrodynamic equations, which describe 
the behavior at large length and time scales, requires a more careful analysis 
of the corresponding Liouville operator jC. Before describing this approach in 
more detail, we provide a more heuristic discussion of the equation of state 
and of one of the transport coefficients, the shear viscosity, using more familiar 
approaches for analyzing the behavior of dynamical systems. 

3.1 

Equation of State 

In a homogeneous fluid, the pressure is the normal force exerted by the fluid 
on one side of a unit area on the fluid on the other side; expressed somewhat 
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differently, it is the momentum transfer per unit axea per unit time across an 
imaginary (flat) fixed surface. There are both kinetic and virial contributions 
to the pressure. The first arises from the momentum transported across the 
surface by particles that cross the surface in the unit time interval; it yields the 
ideal-gas contribution, Pid = NksT/V, to the pressure. For classical particles 
interacting via pair-additive, central forces, the intcrmolccular "potential" 
contribution to the pressure can be determined using the method introduced 
by Irving and Kirkwood [43] . A clear discussion of this approach is given by 
Davis in Rcf. [44], where it is shown to lead to the virial equation of state of 
a homogeneous fiuid. 



in three dimensions, where Fj is the force on particle i due to all the other 
particles, and the sum runs over all particles of the system. 

The kinetic contribution to the pressure, Pid = NksT/V, is clearly present 
in all MFC algorithms. For SRD, this is the only contribution. The reason is 
that the stoc;hastic rotations, which define the c;ollisions, transport (on aver- 
age) no net momentum across a fixed dividing surface. More general MFC 
algorithms (such as those discussed in Sec. 6) have an additional contribution 
to the virial equation of state. However, instead of an explicit force F^ as 
in Eq. (9), the contribution from the multi-particle collisions is a force of the 
form mAvi/At, and the role of the particle position, rj, is played by a variable 
which denotes the cell-partners which participate in the collision [33,45]. 

3.2 

Shear Viscosity 

Just as for the pressure, there are both kinetic and coUisional contributions to 
the transport coeSicients. We present here a heuristic discussion of these con- 
tributions to the shear viscosity, since it illustrates rather clearly the essential 
physics and provides background for subsequent technical discussions. 

Consider a reference plane (a line in two dimensions) with normal in the 
y-direction embedded in a homogeneous fluid in equilibrium. The fluid below 
the plane exerts a mean force Py per unit area on the fluid above the plane; 
by Newton's third law, the fluid above the plane must exert a mean force 
— Pj, on the fluid below the plane. The normal force per unit area is just 
the mean pressure, P, so that Pyy = P. In a homogeneous simple fluid in 
which there are no velocity gradients, there is no tangential force, so that, for 
example, Pyx = 0. Pa/3 is called the pressure tensor, and the last result is just 
a statement of the well known fact that the pressure tensor in a homogeneous 
simple fluid at equilibrium has no off-diagonal elements; the diagonal elements 
are all equal to the mean pressure P. 

Consider a shear flow with a shear rate 7 = dux{y) /dy. In this case, there 
is a tangential stress on the reference surface because of the velocity gradient 




(9) 
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normal to the plane. In the small gradient limit, the dynamic viscosity, r], is 
defined as the coefficient of proportionality between the tangential stress, Py^, 
and the normal gradient of the imposed velocity gradient. 

Pyx = -Vi- (10) 

The kinematic viscosity, v, is related to rj hy v = rj/ p, where p = nm, is the 
mass density, with n the number density of the fluid and m the particle mass. 

Kinetic contribution to the shear viscosity: The kinetic contribution to the 
shear viscosity comes from transverse momentum transport by the flow of 
fluid particles. This is the dominant contribution to the viscosity of gases. 
The following analogy may make this origin of viscosity clearer. Consider two 
ships moving side by side in parallel, but with different speeds. If the sailors 
on the two ships constantly throw sand bags from their ship onto the other, 
there will be a transfer of momentum between to two ships so that the slower 
ship accelerates and the faster ship decelerates. This can be interpreted as an 
effective friction, or kinetic viscosity, between the ships. There are no direct 
forces between the ships, and the transverse momentum transfer originates 
solely from throwing sandbags from one ship to the other. 

A standard result from kinetic theory is that the kinetic contribution to 
the shear viscosity in simple gases is [46] 

7?'=*" ~ nmvX, (11) 

where A is the mean free path and v is the thermal velocity. Using the fact 
that A ~ vAt for SRD and that v <^ ^Jk^TJm, relation (11) implies that 

Ty'^'" ~ nksTAt, or equivalently, ~ kBTAt/m, (12) 

which is, as more detailed calculations presented later will show, the correct 
dependence on n, ksT, and At. In fact, the general form for the kinetic 
contribution to the kinematic viscosity is 

u'''"=''^^fkin{d,M,a), (13) 
m 

where d is the spatial dimension, M is the mean number of particles per 
cell, and a is the SRD collision angle. Another way of obtaining this result 
is to use the analogy with a random walk: The kinematic viscosity is the 
diffusion coefficient for momentum diffusion. At large mean free path. A/a ^ 
1, momentum is primarily transported by particle translation (as in the ship 
analogy). The mean distance a particle streams during one time step. At, 
is A. According to the theory of random walks, the corresponding diffusion 
coefficient scales as i/*"*" ~ X'^/At ~ kBTAt/m. 

Note that in contrast to a "real" gas, for which the viscosity has a square 
root dependence on the temperature, i^*^*" ^ T for SRD. This is because the 
mean free path of a particle in SRD does not depend on density; SRD allows 
particles to stream right through each other between collisions. Note, however, 
that SRD can be easily modified to give whatever temperature dependence 
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is desired. For example, an additional temperature-dependent collision prob- 
ability can be introduced; this would be of interest, e.g., for a simulation of 
realistic shock- wave profiles. 

Collisional contribution to the shear viscosity: At small mean free paths. 
A/a <C 1, particles "stream" only a short distance between collisions, and the 
multi-particle "collisions" are the primary mechanism for momentum trans- 
port. These collisions redistribute momenta within cells of linear size a. This 
means that momentum "hops" an average distance a in one time step, leading 
to a momentum diffusion coefficient I At. The general form of the 

collisional contribution to the shear viscosity is therefore 

^'''' = Xjcoi{d,M,a). (14) 

This is indeed the scaling observed in numerical simulations at small mean 
free path. 

The kinetic contribution dominates for A ^ a, while the collisional con- 
tribution dominates in the opposite limit. Two other transport coefficients of 
interest are the thermal diffusivity, Dt, and the single particle diffusion co- 
efficient, D. Both have the dimension m^/sec. As dimensional analysis would 
suggest, the kinetic and collisional contributions to Dt exhibit the same char- 
acteristic dependencies on A, a, and At described by Eqs. (13) and (14). Since 
there is no collisional contribution to the diffusion coefficient, D ~ X^/At. 

Two complementary approaches have been used to derive the transport co- 
efficients of the SRD fluid. The first is an equilibrium approach which utilizes 
a discrete projection operator formalism to obtain Green-Kubo (GK) relations 
which express the transport coefficients as sums over the autocorrelation func- 
tions of reduced fluxes. This approach was first utilized by Malevanets and 
Kapral [19], and later extended by Ihle, Kroll and Tiizel [20,27,28] to in- 
clude colUsional contributions and arbitrary rotation angles. This approach is 
described in Sec. 4.1. 

The other approach uses kinetic theory to calculate the transport coef- 
ficients in stationary non-equilibrium situation such as shear flow. The first 
application of this approach to SRD was presented in Ref. [21], where the col- 
lisional contribution to the shear viscosity for large M, where particle number 
fluctuations c;an be ignored, was calculated. This scheme was later extended 
by Kikuchi et al. [26] to include fluctuations in the number of particles per 
cell, and then used to obtain expressions for the kinetic contributions to shear 
viscosity and thermal conductivity [35]. This non-equilibrium approach is de- 
scribed in Sec. 5. 

4 

Equilibrium Calculation of Dynamic Properties 

A projection operator formalism for deriving the linearized hydrodynamic 
equations and Green-Kubo (GK) relations for the transport coefficients of 
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molecular fluids was originally introduced by Zwanzig [47-49] and later 
adapted for lattice gases by Dufty and Ernst [50]. With the help of this formal- 
ism, explicit expressions for both the reversible (Euler) as well as dissipative 
terms of the long-time, large-length-scale hydrodynamics equations for the 
coarse-grained hydrodynamic variables were derived. In addition, the result- 
ing GK relations enable explicit calculations of the transport coefficients of 
the fluid. This work is summarized in Sec. 4.1. An analysis of the equilibrium 
fluctuations of the hydrodynamic modes can then be used to directly measure 
the shear and bulk viscosities as well as the thermal diffusivity. This approach 
is described in Sec. 4.2, where SRD results for the dynamic structure factor 
are discussed. 

4.1 

Linearized Hydrodynamics and Green-Kubo Relations 

The Green-Kubo (GK) relations for SRD differ from the well-known continu- 
ous versions due to the discrete-time dynamics, the underlying lattice struc- 
ture, and the multi-particle interactions. In the following, we briefly outline 
this approach for determining the transport coefficients. More details can be 
found in Refs. [20,27]. 

The starting point of this theory are microscopic definitions of local hy- 
drodynamic densities A/}. These "slow" variables are the local number, mo- 
mentum, and energy density. At the cell level, they are defined as 

N d 

A/3(c)=E«wn^(i- 

i=\ 7=1 

with the discrete cell coordinates ^ = am, where mp = 1,...,L, for each 
spatial component, ai^i = 1 is the particle density, {afj^i} = m{wi(^_i)}, with 
/? = 2, d+ 1, are the components of the particle momenta, and ad+2,i = 
mvf/2 is the kinetic energy of particle i. d is the spatial dimension, and rj 
and Vi are position and velocity of particle i, respectively. 

Af-){$^), for [3 = 2, . . . ,d + 2, are cell level coarse-grained densities. For 
example, ^2(^) is the x-component of the total momentuni of all the particles 
in cell ^ at the given time. Note that the particle density, Ai , was not coarse- 
grained in Ref. [20], i.e., the & functions in Eq. (15) were replaced by a 
^-function. This was motivated by the fact that during collisions the particle 
number is trivially conserved in areas of arbitrary size, whereas energy and 
momentum are only conserved at the cell level. 

The equilibrium correlation func;tions for the conserved variables are de- 
fined by (5A^(r, t)5A^{T' , t')), where {6 A) = A— {A), and the brackets denote 
an average over the equilibrium distribution. In a stationary, translationally 
invariant system, the correlation fuiic;tioiis depend only on the differences r— r' 
and t — t', and the Fourier transform of the matrix of correlation functions is 

Ga0{Kt) = :^(M^(k,0)M^(k,t)), (16) 



(15) 
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where the asterisk denotes the complex conjugate, and the spatial Fourier 
transforms of the densities are given by 



A0{\.) = Y,a0,^'^-'^, (17) 
i 

where is the coordinate of the cell occupied by particle j. k = 27rn/(aL) is 
the wave vector, where n/s = 0, ±1, . . . , ±(L— 1), L for the spatial components. 
To simplify notation, we omit the wave-vector dependence of Gai3 in this 
section. 

The collision invariants for the conserved densities are 

y^-e,it+At)^^^.^^ + At) - a0j{t)] = 0, (18) 



where is the coordinate of the cell occupied by particle j in the shifted 
system. Starting from these conservation laws, a projection operator can be 
constructed that projects the full SRD dynamics onto the conserved fields [20]. 
The central result is that the discrete Laplace transform of the linearized 
hydrodynamic equations can be written as 

[s + ikn + A;2yl]G(k, s) = ^G(0)i?(fc), (19) 

where R{k) = [1 + At{ikn + k'^A]~^ is the residue of the hydrodynamic 
pole [20] . The linearized hydrodynamic equations describe the long-time large- 
length-scale dynamics of the system, and are valid in the limits of small k 
and s. The frequency matrix J? contains the reversible (Euler) terms of the 
hydrodynamic equations. A is the matrix of transport coefficients. The discrete 
Green-Kubo relation for the matrix of viscous transport coefficients is [20] 

At °° ' 

= J^jT^ E Ckx(TaXmkx'<T(3X'{t)), (20) 

^ t=0 

where the prime on the sum indicates that the t = term has the relative 
weight 1/2. aai3 = P5ai3—Pa0 is the viscous stress tensor. The reduced fluxes 
in Eq. (20) are given by 

kx^axit) = 5 E (-^.-^ • + ^v^a{t)AeM + fLv^it)) (21) 

j 

for a = 1, . . . , d, with A^j (t) = 1^ {t + At) - (t) , A^j {t + At) = {t + At) - 
^*(t + At), and Avxj{t) = Vxj{t + At) — Vxj{t)- ^j{t) is the cell coordinate 
of particle j at time t, while ^* is it's cell coordinate in the (stochastically) 
shifted frame. The corresponding expressions for the thermal diffusivity and 
self-diffusion coefficient can be found in Ref. [20]. 

The straightforward evaluation of the GK relations for the viscous (21) and 
thermal transport coefficients leads to three — kinetic, collisional, and mixed — 
contributions. In addition, it was found that for mean free paths A smaller 
than the cell size a, there are finite cell-size corrections which could not be 
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summed in a controlled fashion. The origin of the problem was the explicit 
appearance of in the stress correlations. However, it was subsequently 
shown [28,51] that the Grccn-Kubo relations can be rc-summcd by introducing 
a stochastic variable, Bi, which is the difference between change in the shifted 
cell coordinates of particle i during one streaming step and the actual distance 
traveled, AtVi. The resulting microscopic stress tensor for the viscous modes 
is 

= ^ ["I'yiaVi/S + ^ViaBiis (22) 
i 

where Bjp{t) — £,j/j{t + At) — $|^(t) — AtVji}{t). It is interesting to compare 
this result to the corresponding expression 



= ^(5(r - Ti) 



(23) 



for molecular fluids. The first term in both expressions, the ideal-gas contri- 
bution, is the same in both cases. The collisional contributions, however, are 
quit different. The primary reason is that in SRD, the collisional contribution 
corresponds to a nonlocal (on the scale of the cell size) force which acts only 
at discrete time intervals. 

Bi has a number of important properties which simplify the calculation 
of the transport cociEcicnts. In particular, it is shown in Rcfs. [28, 51] that 
stress-stress correlation functions involving one Bi in the GK relations for 
the transport coefficients are zero, so that, for example, A^fi^ = yl^™(k) -|- 

yl-^(k), with 

At °° ' 

'^''^^^^^ = NmkT ^ {h'yltmh"yf;:{nAt)) (24) 

n=0 



and 



At °° 

^^"^^""^ = iVmFr ^ (fcA<°i(0)I^A'<(n^i)>], (25) 
with 

a^^{nAt) = ^mvja{nAt)vji3{nAt) (26) 
j 

and 

<0{n^t) = ^ E mvj^{nAt)Bj0{nAt), (27) 

i 

where Bj0{nAt) = £,ji}{[n + l]At) — ^|^(nZ\f) — AtVjp{nAt). Similar relations 
were obtained for the thermal diffusivity in Ref. [28]. 
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4.1.1 

Explicit Expressions for the Transport Coefficients 

Analytical calculations of the SRD transport coefficients are greatly simplified 
by the fac;t that coUisional and kinetic contributions to the stress-stress auto- 
correlation functions decouple. Both the kinetic and coUisional contributions 
have been calculated expHcitly in two and three dimension, and numerous nu- 
merical tests have shown that the resulting expressions for all the transport co- 
efficients are in excellent agreement with simulation data. Before summarizing 
the results of this work, it is important to emphasize that because of the cell 
structure introduced to define coarse-grained collisions, angular momentum 
is not conserved in a collision [28,35,39]. As a consequence, the macroscopic 
viscous stress tensor is not, in general, a symmetric function of the derivatives 
daVf). Although the kinetic contributions to the transport coefficients lead to 
a symmetric stress tensor, the coUisional do not. Before evaluating the trans- 
port coefficients, we discuss the general form of the macroscopic viscous stress 
tensor. 

Assuming only cubic symmetry and allowing for a non-symmetric stress 
tensor, the most general form of the linearized Navier-Stokes equation is 

5tt;,(k) = -d^p + A^p{\^)vp{\^), (28) 

where 

yl„^(k) = vi ^(5«,/3 + '^-J-f'c.kis^ (29) 

In a normal simple liquid, k — (because of invariance with respect to in- 
finitesimal rotations) and 1^2 = (because the stress tensor is symmetric in 
daVp), SO that the kinematic shear viscosity is u = vi. In this case, Eq. (29) 
reduces to the well-known form [20] 

Aafsiii) = V [5a,f3 + '^-^kakff^ + -fkak/B , (30) 

where 7 is the bulk viscosity. 

Kinetic contributions: Kinetic contributions to the transport coefficients dom- 
inate when the mean free path is larger than the cell size, i.e., X > a. As can 
be seen from Eqs. (24) and (26), an analytic calculations of these contribu- 
tions requires the evaluation of time correlation functions of products of the 
particle velocities. This is straightforward if one makes the basic assumption 
of molecular chaos that successive collisions between particles are not corre- 
lated. In this case, the resulting time-series in Eq. (24) is geometrical, and 
can be summed analytically. The resulting expression for the shear viscosity 
in two dimensions is 



2m V(M-l + e-^)sin2(a) ^ ' 
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Fig. 1. (a) Normalized kinetic contribution to the viscosity, u''^" / (AtkBT), in three 
dimensions as a function of the collision angle a. Data were obtained by time aver- 
aging the Green-Kubo relation over 75,000 iterations using A/a = 2.309 for M = 5 
(■) and M = 20 (•). The lines are the theoretical prediction, Eq. (32). Parameters: 

L/a = 32, Z\t = 1. From Ref. [53]. 

(b) Normalized coUisional contribution to the viscosity, v'^°'' At / , in three dimen- 
sions as a function of the collision angle a. The solid line is the theoretical predic- 
tion, Eq. (39). Data were obtained by time averaging the Green-Kubo relation over 
300,000 iterations. Parameters: L/a = 16, X/a = 0.1, M = 3, and At = 1. Prom 
Ref. [54]. 



Fluctuations in the number of particles per cell are included in (31). This 
result agrees with the non-equilibrium calculations of Pooley et al. [35,52], 
measurements in shear flow [26], and the numerical evaluation of the GK 
relation in equilibrium simulations (see Fig. 1). 

The corresponding result in three dimensions for collision rule (3) is 

Mn ksTAt { 5M 



2m V (M - 1 + e-^)[2 - cos(a) - cos(2a)] 7 " ^^^^ 

The kinetic contribution to the stress tensor is symmetric, so that z/f'" = 
and the kinetic contribution to the shear viscosity is z^*^*" = z/f*". 
CoUisional contributions: Explicit expressions for the colUsional contributions 
to the viscous transport coefficients can be obtained by considering various 
choices for k and a and (3 in Eqs. (25), (27) and (29). Taking k in the y- 
direction and a = (5 = 1 yields 



t=0 i,3 

Other choices lead to relations between the coUisional contributions to the 
viscous transport coefficients, namely 

[l + {d-2)/<]\ ul"^ + 7™' + k""^ = vf^ + . (34) 

and 

[{d - 2) /d]vl°^ - vl°^ + 7^"' = 0. (35) 
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These results imply that k™' = 0, and 7™' - 2i/f '/d = - It follows 
that the collision contribution to the macroscopic viscous stress tensor is 

^lilP = <{dpv^ + d^vp) + vf{d0Vo, - davp) + [uf - yT')5„0dxvx 

= « + U^"')dfSVc, + {V^°' - (36) 

where Qajs = Sa0d\v\ — daVp. Since Qajs has zero divergence, djsQa^ = 0, the 
term containing Q in Eq. (36) will not appear in the linearized hydrodynamic 
equation for the momentum density, so that 

= -Vp + p(!/'=^" + v''°^)Aw + ^^!/'=^"V(V • v), (37) 

where v'^°^ = vf^ + ff°'. In writing Eq. (37) we have used the fact that the 
kinetic contribution to the microscopic stress tensor, ct*^™, is symmetric, and 
^km _ Q J27] . The viscous contribution to the sound attenuation coefficient is 
yColj^2(d- l)i/''*"/fi instead of the standard result, 2{d-l)v/d+^, for simple 
isotropic fluids. The coUisional contribution to the effective shear viscosity is 
ycol = ycoi _|_ j^coj is interesting to note that the kinetic theory approach 
discussed in Ref. [35] is able to show explicitly that z/™' = 1/2°') so that 

It is straightforward to evaluate the various contributions to the right hand 
side of (33). In particular, note that since velocity correlation functions are 
only required at equal times and for a time lag of one time step, molecular 
chaos can be assumed [51] . Using the relation [28] 

{Bi„{nAt)Bji3{mAt)) = — (5^/3(1 + 5^) [2(5„,„ - 5n,m+i - ^n,m-i] , (38) 

and averaging over the number of particles in a cell assuming that the number 
of particles in any cell is Poisson distributed at each time step, with an average 
number M of particles per cell, one then finds 

for the SRD collision rules in both two and three dimensions. Eq. (39) agrees 
with the result of Refs. [26] and [35] obtained using a completely different 
non-equilibrium approach in shear flow. Simulation results for the collisional 
contribution to the viscosity are in excellent agreement with this result (see 
Fig. 1). 

Thermal diffusivity and self-diffusion coefficient: As with the viscosity, there 
are both kinetic and collisional contributions to the thermal diffusivity, Dt- 
A detailed analysis of both contributions is given in Ref. [28] , and the results 
are summarized in Table 1. The self-diffusion coefficient, D, of particle i is 
defined by 

1 At °° ' 

n=0 
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where the second expressions is the corresponding discrete GK relation. The 
self-diffusion coefficient is unique in that the collisions do not explicitly con- 
tribute to D. With the assumption of molecular chaos, the kinetic contribu- 
tions are easily summed [27] to obtain the result given in Table 1. 

4.1.2 

Beyond Molecular Chaos 

The kinetic contributions to the transport coefficients presented in Table 1 
have all been derived under the assumption of molecular chaos, i.e., that 
particle velocities are not correlated. Simulation results for the shear viscosity 

and tlK^rnial diffusivity have generally been found to be in good agreement 
with these results. However, it is known that there are correlation effects for 
A/a smaller than unity [15,55]. They arise from correlated collisions between 
particles that arc in the same collision cell for more than one time step. 

For the viscosity and thermal conductivity, these corrections are generally 
negligible, since they are only significant in the small A/a regime, where the 
coUisional c;ontribution to the transport coefficients dominates. In this regard, 
it is important to note that there are no correlation corrections to 1/"^°' and 
I?™' [28]. For the self-diffusion coefficient — for which there is no collisional 
contribution correlation corrections dramatically increase the value of this 
transport coefficient for A <C a, see Refs. [15,55]. These correlation corrections, 
which arise from particles which collide with the same particles in consecutive 
time steps, are distinct from the correlations effects which are responsible for 
the long-time tails. This distinction is important, since long-time tails are also 
visible at large mean free paths, where these corrections are negligible. 

4.2 

Dynamic Structure Factor 

Spontaneous thermal ffuctuations of the density, p{r,t), the momentum den- 
sity, g{r,t), and the energy density, e{r,t), are dynamically coupled, and an 
analysis of their dynamic correlations in the limit of small wave numbers and 
frequencies can be used to measure a fluid's transport coefficients. In par- 
ticular, because it is easily measured in dynamic light scattering, x-ray, and 
neutron scattering experiments, the Fourier transform of the density-density 
correlation function — the dynamics structure factor — is one of the most widely 
used vehicles for probing the dynamic and transport properties of hquids [56] . 

A detailed analysis of equilibrium dynamic correlation functions the dy- 
namic structure factor as well as the vorticity and entropy-density correlation 
functions — using the SRD algorithm is presented in Ref. [57]. The results — 
which are in good agreement with earlier numerical measurements and the- 
oretical predictions — provided further evidence that the analytic expressions 
or the transport coefficients are accurate and that we have an excellent un- 
derstanding of the SRD algorithm at the kinetic level. 
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d 


Kinetic {xkBTAt/2m) 


Collisional 


(xa^/Zit) 




2 
3 


M 1 
(M-l + e-M) sin2(Q,) 

5M -1 

(M-l+e-*')[2-coB(a)-cos(2a)] 


(M-l + e-*f) 
6dM 


[1 — cos(a)] 


Dt 


2 
3 




(1-1/M) r, 
3(d+2)M L-^ 


— cos(a)] 


D 


2 
3 


dM 1 

(l-cos(a))(M-l+e-'^) 





Table 1. Theoretical expressions for the kinematic shear viscosity p, the thermal 
diffusivity, Dt, and the self-difTusion coefficient, D, in both two (d = 2) and three 
(d = 3) dimensions. M is the average number of particles per cell, a is the collision 
angle, fcs is Boltzmann's constant, T is the temperature, At is the time step, m is 
the particle mass, and a is the cell size. Except for self-diffusion constant, for which 
there is no collisional contribution, both the kinetic and collisional contributions are 
listed. The expressions for shear viscosity and self-diffusion coefficient include the 
effect of ffuctuations in the number of particles per cell; however, for brevity, the 
relations for thermal diffusivity are correct only up to 0{1/M) and 0(1/M^) for the 
kinetic and collisional contributions, respectively. For the complete expressions, see 
Refs. [28,53,54]. 



Here, we briefly summarize the results for the dynamic structure factor. 
The dynamic structure exhibits three peaks, a central "Rayleigh" peak caused 
by the thermal diffusion, and two symmetrically placed "Brillouin peaks" 
caused by sound. The width of the central peak is determined by the thermal 
diffusivity, Dt, while that of the two Brillouin peaks is related to the sound 
attenuation coefficient, F. For the SRD algorithm [57], 

r = DT(^^-l) +2(^^^z/'=™ + i/^°'. (41) 

Note that in two-dimensions, the sound attenuation coefficient for a SRD fluid 
has the same functional dependence on Dt and f = v'^^^ + 1/"°'- as an isotropic 
fluid with an ideal-gas equation of state (for which 7 = 0). 

Simulation results for the structure factor in two-dimensions with A/a = 
1.0 and collision angle a = 120°, and A/a = 0.1 with collision angle a = 60° 
are shown in Figs. 2a and 2b, respectively. The solid Unes are the theoretical 
prediction for the dynamic structure factor (sec Eq. (36) of Rcf. [57]) using 
c = y/2kBT/m and values for the transport coefficients obtained using the 
expressions in Table 1, assuming that the bulk viscosity 7 = 0. As can be 
seen, the agreement is excellent. 
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Fig. 2. Normalized dynamic structure, S'pp{kLo) / Xpp{k) , for k = 27r(l, and (a) 
A/a = 1.0 with a = 120°, and (b) A/a = 0.1 with a = 60°. The sohd hncs are the 
theoretical prediction for the dynamic structure factor (see Eq. (36) of Ref. [57]) 
using values for the transport coefficients obtained using the expressions in Table 1. 
The dotted lines show the predicted positions of the Brillouin peaJcs, w = ±cfc, with 
c = -^2fesT/m. Parameters: L/a = 32, M = 15, and At = 1.0. Prom Ref. [57]. 



5 

Non-Equilibrium Calculations of Transport Coefficients 

MPC transport coefficients have also be evaluated by calculating the linear 
response of the system to imposed gradients. This approach was introduced 
by Kikuchi et al. [26] for the shear viscosity and then extended and refined 
in Ref. [35] to determine the thermal diffusivity and bulk viscosity. Here, we 
summarize the derivation of the shear viscosity. 



5.1 

Shear Viscosity of SRD: Kinetic Contribution 

Linear response theory provides an alternative, and complementary, approach 
for evaluating the shear viscosity. This non-equilibrium approach is related 
to equilibrium calculations described in the previous section through the 
fluctuation-dissipation theorem. Both methods yield identical results. For the 
more complicated analysis of the hydrodynamic equations, the stress tensor, 
and the longitudinal transport coefficients such as the thermal conductivity, 
the reader is referred to Ref. [35]. 

Following Kikuchi et al. [26], we consider a two-dimensional liquid with 
an imposed shear 7 = dux{y)/dy. On average, the velocity profile is given 
by V = (7;/;, 0). The dynamic shear viscosity rj is the proportionality constant 
between the velocity gradient 7 and the frictional force acting on a plane 
perpendicular to y; i.e. 

<^xy = f?7, (42) 

where cjxy is the off-diagonal element of the viscous stress tensor. During the 
streaming step, particles will cross this plane only if \vy At\ is greater than the 
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distance to the plane. Assuming that the fluid particles are homogeneously 
distributed, the momentum flux is obtained by integrating over the coordi- 
nates and velocities of all particles that cross the plane from above and below 
during the time step At. The result is [26] 



'xy 



7^* /..2x 



«) - M , (43) 



where the mass density p = mM/a'^, and the averages are taken over the 
steady-state distribution P{vx — 7y,%). It is important to note that this is 
not the Maxwell-Boltzmann distribution, since we are in a non-equilibrium 
steady state where the shear has induced correlations between and Vy. As 

a consequence, (vxVy) is nonzero. To determine the behavior of [v^Vy), the 
effect of streaming and collisions are calculated separately. During streaming, 
particles which arrive at yo with positive velocity Vy have started from yo — 

Vy At; these particles bring a velocity component which is smaller than that 
of particles originally located at yo. On the other hand, particles starting out at 
y > yo with negative Vy bring a larger Vx ■ The velocity distribution is therefore 
sheared by the streaming, so that P"^*'^^{vx,Vy) = P^^^°'^^[vx + 'yVyAt,Vy). 
Averaging v^Vy over this distribution gives [26] 

i^xvyr^"'^ = (v^^vy) - jAti^^ , (44) 

where the superscript denotes the quantity after streaming. The streaming 
step therefore reduces correlations by —'yAt{Vy), making Vx and Vy increas- 
ingly anti-corrclatcd. 

The collision step redistributes momentum between particles and tends to 
reduce correlations. Making the assumption of molecular chaos, i.e., that is 
that the velocities of different particles arc uncorrelated, and averaging over 
the two possible rotation directions, one finds. 



{VxVy) 



after 



l-^V^[l-cos(2a)] 



{VxVy)'''^'"'' (45) 



The number of particles in a cell, Nc is not constant, and density fluctuations 
have to be included. The probability to find n uncorrelated particles in a 
given cell is given by the Poisson distribution, w{n) = exp(— M)M"/n!; the 
probability of a given particle being in a cell together with n — 1 others is 
nw{n)/M. Taking an average over this distribution gives 

{VxVyTf"^' = f{VxVy)'''f''^'^, (46) 

with 



f 



M - 1 + cxp(-M) 



[1 - cos(2a)] 



(47) 



The difference between this result and just replacing Nc by M in Eq. (45) is 
small, and only important for M < 3. One sees that (vxVy) is first modified by 
streaming and then multiplied by a factor / in the subsequent collision step. 
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In the steady state, it therefore oscillates between two values. Using Eqs. (44), 
(46), and (47), we obtain the self-consistency condition {{v^Vy) —^At{Vy))f = 
{vxVy). Solving for {v^Vy), assuming equipartition of energy, {Vy) = kBT/m, 
and substituting into (43), wc have 

jMAtknTfl f \ 
^^v = ( o + 1 ) ' (48) 



m \2 1-f 

Inserting this result into the definition of the viscosity, (42), yields the same 

expression for the kinetic viscosity in two-dimensional as obtained by the 
equilibrium Green-Kubo approach discussed in Sec. 4.1.1. 



5.2 

Shear Viscosity of SRD: Collisional Contribution 

The collisional contribution to the shear viscosity is proportional to a'^ /At; as 
discussed in Sec. 3.2, it results from the momentum transfer between particles 
in a cell of size a during the colhsion step. Consider again a collision cell of 
linear dimension a with a shear flow Ux{y) = '^/y. Since the collisions occur in 
a shifted grid, they cause a transfer of momentum between neighboring cells 
of the original unshifted reference frame [21,27]. Consider now the momen- 
tum transfer due to collisions across the line y = h, the coordinate of a cell 
boundary in the unshifted frame. If we assume a homogeneous distribution 
of particles in the collision cell, the mean velocities in the upper {y > h) and 
lower partitions are 



"i = ]i^E^^ U2 = ^ 5^ V,, (49) 

^ j=l ^ i=Mi + l 

respectively, where Mi = M{a — h)/a and M2 = Mh/a. Collisions trans- 
fer momentum between the two parts of the cell. The a;-component of the 

momentum transfer is 

Ml 

Ap,{h)^^[vtt''-vt^n . (50) 

i=l 

The use of the rotation rule (2) together with an average over the sign of the 
stochastic rotation angle yields 

^Px(/i) = [cos(a)-l]Mi(ui^-w^) . (51) 

Since Mu = Miui -|- M2U2, 

Apx{h) = [1 - cos(a)] M (u^x - uix) - ( 1 - -] . (52) 

a \ a J 

Averaging over the position h of the dividing line, which corresponds to av- 
eraging over the random shift, we find 
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If 1 

{Ap^) = - / Ap^{h)dh = ^[1 - cos{a)]M{u2x - ui^)- (53) 
d Jo " 



Since the dynamic viscosity r] is defined as the ratio of the tangential stress, 
Pyx, to dux/dy, we have 

dux/dy {u2x - uix)/{a/2)' 

so that the kinematic viscosity, v = rj/ p, in two-dimensions for SRD is 

v'''' = ^^{l-cos{a)] (55) 

in the limit of small mean free path. Since we have neglected the fluctuations 
in the particle number, this expression corresponds to the limit M ^ oo. Even 
though this derivation is somewhat heuristic, it gives a remarkably accurate 

expression; in particular, it contains the correct dependence on the cell size, 
a, and the time step. At, in the limit of small free path, 

^'"' = Xjcoi{d,M,a), (56) 

as expected from simple random walk arguments. Kikuchi et al. [26] included 
particle number fluctuations and obtained id(mtical results for the coUisional 
contribution to the viscosity as was obtained in the Green-Kubo approach 
(see Table 1). 



5.3 

Shear Viscosity of iVIPC-AT 

For MPC-AT, the viscosities have been calculated in Ref. [32] using the meth- 
ods described in Sees. 5.1 and 5.2. The total viscosity of MPC-AT is given by 

the sum of two terms, the coUisional and kinetic contributions. For MPC-AT- 
a, it was found for both two and three dimensions that [32] 

ksTAt [ M 1 



^'-.2aA M )■ 

The exponential terms are due to the fluctuation of the particle num- 
ber per cell and become important for M < 3. As was the case for SRD, 
the kinetic viscosity has no anti-symmetric component; the coUisional con- 
tribution, however, does. Again, as discussed in Sec. 4.1.1 for SRD, one finds 
z/^oi _ jycoi _ ycoi j2. This relation is true for all —a versions of MPC discussed 
in Refs. [32,58,59]. Simulation results were found to be in good agreement 
with theory. 

For MPC- AT -Fa it was found for sufliciently large M that [38, 59] 
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ksTAt ( M 1 



TO \M -{d+lL)l^ 2 



, ^ / M-7/5 \ 
"24Zitlv j- ^^^^ 

MPC-AT-a and MPC-AT+a both have the same kinetic contribution to the 
viscosity in two dimensions; however, imposing angular-momentum conserva- 
tion makes the collisional contribution to the stress tensor symmetric, so that 
the asymmetric contribution, v^, discussed in Sec. 4.1.1 vanishes. The result- 
ing collisional contribution to the viscosity is then reduced by a factor close 
to two. 



6 

Generalized MPC Algorithms for Dense Liquids and Binary 
Mixtures 

The original SRD algorithm models a single-component fluid with an ideal- 
gas equation of state. The fluid is therefore very compressible, and the speed 
of sound, Cs, is low. In order to have negligible compressibility effects, as in 
real liquids, the Mach number has to be kept small, which means that there 
are limits on the flow velocity in the simulation. The SRD algorithm can be 
modified to model both excluded volume effects, allowing for a more realistic 
modeling of dense gases and liquids, as well as repulsive hard-core interac- 
tions between components in mixtures, which allow for a thermodynamically 
consistent modeling of phase separating mixtures. 



6.1 

Non-Ideal Model 

As in SRD, the algorithm consists of individual streaming and collision steps. 
In order to define the collisions, a second grid with sides of length 2a is intro- 
duced, which (in d = 2) groups four adjacent cells into one "supercell". The 
cell structure is sketched in Fig. 3 (left panel). To initiate a collision, pairs 
of cells in every supercell are chosen at random. Three different choices are 
possible: a) horizontal (with <ti = x), b) vertical [a^ = y), and c) diagonal 
collisions (with cr^ = {x + y)/\/2 and (T4 = (x — y)/\/2). 

For a mean particle velocity u„ = (1/M„) of cell n, the projec- 

tion of the difference of the mean velocities of the selected cell pairs on aj, 
Au = (Tj ■ (ui — U2), is then used to determine the probability of collision. If 
Au < 0, no collision will be performed. For positive Au, a collision will occur 
with an acceptance probability, pa, which depends on Au and the number 
of particles in the two cells, Ni and 7V2. The choice of pA determines both 
the equation of state and the values of the transport coefficients. While there 
is considerable freedom in choosing pA, the requirement of thermodynamic 
consistency imposes certain restrictions [33,34,55]. One possible choice is 
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Fig. 3. (left panel) Schematic of collision rules. Momentum is exchanged in three 
ways: a) horizontally along cri, b) vertically along (T2, and c) diagonally along era 
and (74. w and Wd denote the probabilities of choosing collisions a), b), and c) 
respectively. 

(right panel) Static structure factor ^(fe, t = 0) as a function of At for M = 3. 
The open circles (o) show results obtained by taking the numerical derivative of 
the pressure. The bullets (•) are data obtained from direct measurements of the 
density fluctuations. The solid line is the theoretical prediction obtained using the 
first term in Eq. (61) and Eq. (63). fe is the smallest wave vector, fe = (27r/Z/)(l, 0). 
Parameters: L/a = 32, ^ = 1/60, and fesT = 1.0. From Ref. [33]. 



PA(Mi,M2,/iu) = 6>(/iw) tanh(^) with A = AAuNxN^, (59) 

where Q is the unit step function and A is a parameter which is used to tune 
the equation of state. The choice A ^ 'Ni'N^ leads to a non- ideal contribution 
to the pressure which is quadratic in the particle density. 

The collision rule chosen in Ref. [33] maximizes the momentum trans- 
fer parallel to the connecting vector a j and docs not change the transverse 
momentum. It exchanges the parallel component of the mean velocities of 
the two cells, which is equivalent to a "reflection" of the relative velocities, 
v\(t + At~) — fi" = —(y\(t~) — u"), where u" is the parallel component of the 
mean velocity of the particles of both cells. This rule conserves momentum 
and energy in the cell pairs. 

Because of x — y symmetry, the probabilities for choosing cell pairs in the 
X- and y- directions (with unit vectors ai and (72 in Fig. 3) axe equal, and 
will be denoted by w. The probability for choosing diagonal pairs {as and 
0-4 in Fig. 3) is given by Wd = I — 2w. w and must be chosen so that the 
hydrodynamic equations are isotropic and do not depend on the orientation of 
the underlying grid. An equivalent criterion is to guarantee that the relaxation 
of the velocity distribution is isotropic. These conditions require w = 1/4 and 
Wd = 1/2. This particular choice also ensures that the kinetic part of the 
viscous stress tensor is isotropic [45]. 
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6.1.1 

Transport Coefricients 

The transport coefficients can be determined using the same Green-Kubo for- 
malism as was used for the original SRD algorithm [21,51]. Alternatively, the 
non-equiHbrium approach describe in Sec. 5 can be used. Assuming molecular 
chaos and ignoring fluctuations in the number of particles per cell, the kinetic 
contribution to the viscosity is found to be 



which is in good agreement with simulation data, pcoi is essentially the collision 
rate, and can be obtained by averaging the acceptance probability, Eq. (59). 
The coUisional contribution to the viscosity is = Pco;('^^/3^^) [60]. 

The sclf-difFusion constant, D, is evaluated by summing over the velocity- 
autocorrelation function (see, e.g. Ref. [21]); which yields D = Vhin- 

6.1.2 

Equation of State 

The collision rules conserve the kinetic energy, so that the internal energy 
should be the same as that of an ideal gas. Thermodynamic consistency there- 
fore requires that the non-ideal contribution to the pressure is linear in T. This 
is possible if the coefficient A in Eq. (59) is sufficiently small. 

The mechanical definition of pressure — the average longitudinal momen- 
tum transfer across a fixed interface per unit time and unit surface area — can 
be used to determine the equation of state. Only the momentum transfer 
due to collisions needs to be considered, since that coming from streaming 
constitutes the ideal part of the pressure. Performing this calculation for a 
fixed interface and averaging over the position of the interface, one finds the 
non-ideal part of the pressure. 



Pn is quadratic in the particle density, p = M/a^, as it would be expected 

from a virial expansion. The prcfactor A must be chosen small enough that 
higher order terms in this expansion are negligible. Prefactors A leading to 
acceptance rates of about 15% are sufficiently small to guarantee that the 
pressure is linear in T. 

The total pressure is the average of the diagonal part of the microscopic 
stress tensor. 




with 




(60) 




(61) 




(62) 
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The first term gives the ideal part of the pressure, Pid, as discussed in Ref. [21]. 
The average of the second term is the non- ideal part of the pressure, P„. z*; is a 
vector which indexes collision partners. The first subscript denotes the particle 
number and the second, /, is the index of the collision vectors ai in Fig. 3 (left 
panel). The components of z*, are either 0, 1, or —1 [55]. Simulation results 
for Pn obtained using Eq. (62) arc in good agreement with the analytical 
expression, Eq. (61). In addition, measurements of the static structure factor 
S{k — > 0, f = 0) agree with the thermodynamic prediction 

S{k ^ 0, i = 0) = pkBTdp/dP\T (63) 

when result (61) is used [see Fig. 3 (right panel)]. The adiabatic speed of sound 

obtained from simulations of the dynamic structure factor is also in good 
agreement with the predictions following from Eq. (61). These results provide 
strong evidence for the thermodynamic consistency of the model. Consistency 
checks arc particularly important because the non-ideal algorithm docs not 
conserve phase-space volume. This is because the collision probability depends 
on the difference of collision-cell velocities, so that two different states can be 
mapped onto the same state by a collision. While the dynamics presumably 
still obeys detailed — or at least semi-detailed — balance, this is very hard to 
prove, since it would require knowledge not only of the transition probabilities, 
but also of the probabilities of the individual equilibrium states. Nonetheless, 
no inconsistencies due to the absence of time-reversal invariance or a possible 
violation of detailed balance have been observed. 

The structure of S{k) for this model is also very similar to that of a simple 
dense fluid. In particular, for flxed M, both the depth of the minimum at 
small k and the height of the first peak increase with decreasing At, until 
there is an order- disorder transition. The four-fold symmetry of the resulting 
ordered state — in which clusters of particles are concentrated at sites with the 
periodicity close, but not necessarily equal, to that of the underlying grid — 
is clearly dictated by the structure of the collision cells. Nevertheless, these 
ordered structures are similar to the low-temperature phase of particles with 
a strong repulsion at intermediate distances, but a soft repulsion at short 
distances. The scaling behavior of both the self-diffusion constant and the 
pressure persists until the order/disorder transition. 

6.2 

Phase-Separating Multi-Component Mixtures 

In a binary mixture of A and B particles, phase separation can occur when 
there is an effective repulsion between A-B pairs. In the current model, this is 
achieved by introducing velocity-dependent multi-particle collisions between 
A and B particles. There are Na and Nb partic;les of type A and B, respec- 
tively. In two dimensions, the system is coarse-grained into {L/a)"^ cells of a 
square lattice of linear dimension L and lattice constant a. The generalization 
to three dimensions is straightforward. 
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Collisions are defined in the same way as in the non-ideal model discussed 
in the previous section. Now, however, two types of collisions are possible for 
each pair of cells: particles of type A in the first cell can undergo a collision 
with particles of type B in the second cell; vice versa, particles of type B in 
the first cell can undergo a collision with particles of type A in the second 
cell. There arc no A-A or B-B collisions, so that there is an effective repulsion 
between A-B pairs. The rules and probabilities for these collisions are chosen in 
the same way as in the non-ideal single-component fluid described in Refs. [33, 
55]. For example, consider the collision of A particles in the first cell with the 
B particles in the second. The mean particle velocity of A particles in the 
first cell is u^i — (I/Nc^a) X^il^i^ '^i' where the sum runs over all A particles, 
Nc^A, in the first cell. Similarly, = (1/A^c,b) ^i=f is the mean velocity 
of B particles in the second cell. The projection of the difference of the mean 
velocities of the selected cell-pairs on Uj, Auab = '^j ' (ua ~ u^), is then 
used to determine the probability of collision. If Auab < 0, no collision will 
be performed. For positive Auab, a collision will occur with an acceptance 
probability 



where O is the unit step function and A is a parameter which allows us to 

tunc the equation of state; in order to ensure thermodynamic consistency, it 
must be sufficiently small that that pA < ^ for essentially all collisions. When 
a collision occurs, the parallel component of the mean velocities of colliding 
particles in the two cells, + At) — u\q = —{vf{t) — u\^), is exchanged, 

where u\g = {Nc^aUa + Nc.bu'b) li^c,A + -^c,b) is the parallel component 
of the mean velocity of the colliding particles. The perpendicular component 
remains unchanged. It is easy to verify that these rules conserve momentum 
and energy in the cell pairs. The collision of B particles in the first cell with 
A particles in the second is handled in a similar fashion. 

Because there are no A-A and B-B collisions, additional SRD collisions at 
the cell level are incorporated in order to mix particle momenta. The order of 
A-B and SRD collision is random, i.e., the SRD collision is performed first with 
a probability 1/2. If necessary, the viscosity can be tuned by not performing 
SRD collisions every time step. The results presented here were obtained using 
a SRD collision angle oi a — 90°. 

The transport coefficients can be calculated in the same way as for the 
one-component non-ideal system. The resulting kinetic contribution to the 
viscosity is 



where Ma = (A^cyi), Mb = {Nc.b)- In deep rjucnchcs, the cxnicx^ntration 
of the minority component is very small, and the non-ideal contribution to 
the viscosity approaches zero. In this case, the SRD collisions provide the 
dominant contribution to the viscosity. 



Pa{N^^a.N^^b,Auab) ^ AAuabO{Auab)Nc,aNc,b 



(64) 




(65) 
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6.2.1 

Free Energy 

An analytic expression for the equation of state of this model can be derived 
by calculating the momentum transfer across a fixed surface, in much the 
same way as was done for the non-ideal model in Ref. [33]. Since there are 
only non-ideal collisions between A-B particles, the resulting contribution to 
the pressure is 

Pn=(w+ AMaMb^ = rpAPB, (66) 

where pa and pB are the densities of A and B and F = {10 + Wd/y/2)a^A/At. 
In simulations, the total pressure can be measured by taking the ensemble 
average of the diagonal components of the inic;rosc;opic; stress tensor. In this 
way, the pressure can be measured locally, at the cell level. In particular, the 
pressure in a region consisting of Nceii cells is 

= aJi^, (£ - ^-ix^J^]) , (67) 

-—1 jGc I 



where the second sum runs over the particles in cell c. The first term in 
Eq. (67) is the ideal-gas contribution to the pressure; the second comes from 
the momentum transfer between cells involved in the collision indexed by 
[45]. 

Expression (66) can be used to determine the entropy density, s. The 
ideal-gas contribution to s has the form [61] 

Sideal = PV{T) - kB [pa In PA + Ps In pb] , (68) 

where p = pa + Pb- Since ^p{T) is independent of pA and pB, this term does 
not play a role in the current discussion. The non-ideal contribution to the 
entropy density, s„, can be obtained from Eq. (66) using the thermodynamic 
relation 

Pn/T = -Sn + PAdSn/dpA + psdSn/dpB- (69) 

The result is s„ = Fp^PB, so that the total configurational contribution to 
the entropy density is 

s = -fcs {pAlnpA + PslnpB + Tp^Pb} • (70) 

Since there is no configurational contribution to the internal energy in this 
model, the mean-field phase diagram can be determined by maximizing the 
entropy at fixed density p. The resulting demixing phase diagram as a function 
of Pab = {pa — Pb)/p is given by the sohd line in Fig. 4 (left panel). The 
critical point is located at pab = 0, pF* = 2. For pF < 2. the order parameter 
pAB = 0; for pF > 2, there is phase separation into coexisting A- and B-rich 
phases. As can be seen, the agreement between the mean-field predictions and 
simulation results is very good except close to the critical point, where the 
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Fig. 4. (left panel) Binary phase diagram. There is phase separation for pF > 2. Sim- 
ulation results for pab obtained from concentration histograms are shown as bullets. 
The dashed line is a plot of the leading singular behavior, pab = \/2i{pr — 2)/2, 
of the order parameter at the critical point. The inset shows a configuration 50,000 
time steps after a quench along pab = to pF = 3.62 (arrow). The dark (blue) 
and light (white) spheres are A and B particles, respectively. Parameters: L/a — 64, 
Ma = Mb= 5, fcsT = 0.0004, At = 1, and a = 1. From Ref. [45]. 
(right panel) Dimensionless radial fluctuations, as a function of the mode 

number fe for 4 = 0.45 (•) and A = 0.60 (o) with fcsT = 0.0004. The average 
droplet radii are ro = 11.95 a and ro = 15.21a, respectively. The solid lines are 
fits to Eq. (71). The inset shows a typical droplet configuration for pab = —0.6, 
pF = 3.62 {A = 0.60 and fcsT = 0.0004). Parameters: L/a = 64, Ma = 2, Mb = 8, 
At = 1, a.nd a = 1. From Ref. [45]. 

histogram method of determining the coexisting densities is unrehable and 
critical fluctuations influence the shape of the coexistence curve. 



Surface Tension 

A typical configuration for pab = 0, pP = 3.62 is shown in the inset to Fig. 4 
(left panel), and a snapshot of a fluctuating droplet at pab — —0.6, pF — 3.62 
is shown in the inset to Fig. 4 (right panel). The amplitude of the capillary 
wave fluctuations of a droplet is determined by the surface tension, a. Using 
the parameterization r{(j)) = tq [l -I- J^kL-oo Uk exp{ik(j))'\ and choosing uq to 
fix the area of the droplet, it can be shown that [54] 



Figure 4 (right panel) contains a plot of as a function of mode number 

k for pF = 3.62 and pF = 2.72. Fits to the data yield a ~ 2.9 fcijT for pF = 
3.62 and a ~ 1.1 /c^T for pF = 2.72. Mechanical equilibrium requires that 



6.2.2 




(71) 
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the pressure difference across the interface of a droplet satisfies the Laplace 
equation 

Ap = Pin- Pout ^ id-l)a/ro (72) 

in d spatial dimensions. Measurements of Ap [using Eq. (67)] as a function 
of the droplet radius for A = 0.60 at ksT — 0.0005 yield results in excel- 
lent agreement with the Laplace equation for the correct value of the surface 
tension [45]. 




Fig. 5. (left panel) Droplet configuration in a mixture with Na = 8192, Nb = 32768, 
and Nd — 1500 dimers after 10^ time steps. The initial configuration is a droplet 
with a homogeneous distribution of dimers. The dark (blue) and light (white) colored 
spheres indicate A and B particles, respectively. For clarity, A particles in the bulk 
are smaller and B particles in the bulk are not shown. Parameters: L/a — 64, 
Ma = 2, Mb ^8, A = 1.8, kBT = 0.0001, At = 1, and a = 1. 

(right panel) Typical configuration showing the bicontinuous phase for Na = Nb = 
20480 and Nd = 3000. Parameters: L/a = 64, Ma = 5, Ms = 5, A = 1.8, fcsT = 
0.0001, At = 1, and a = 1. From Ref. [45]. 



The model therefore displays the correct thermodynamic behavior and 
interfacial fluctuations. It can also be extended to model amphiphilic mixtures 
by introducing dimers consisting of tethered A and B particles. If the A and B 
components of the dimers participate in the same collisions as the solvent, they 
behave like amphiphilic molecules in binary oil-water mixtures. The resulting 
model displays a rich phase behavior as a function of pF and the number of 
dimers, Nd- Both the formation of droplets and micelles, as shown in Fig. 5 
(left panel), and a bicontinuous phase, as illustrated in Fig. 5 (right panel), 
have been observed [45] . The coarse-grained nature of the algorithm therefore 
enables the study of large time scales with a feasible computational effort. 
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6.2.3 

Color Models for Immiscible Fluids 

There have been other generaUzations of SRD to model binary mixtures by 
Hashimoto et al. [62] and Inoue et al. [63], in which a color charge, Ci = ±1 
is assigned to two different species of particles. The rotation angle a in the 
SRD rotation step is then chosen such that the color-weighted momentum in 
a cell, m = Y^fl^ Ciiyi — u), is rotated to point in the direction of the gradient 
of the color field c — X]i=i ^i- This rule also leads to phase separation. Sev- 
eral tests of the model have been performed; Laplace's equation was verified 
numerically, and simulation studies of spinodal decomposition and the defor- 
mation of a falling droplet were performed [62]. Later applications include a 
study of the transport of slightly deformed immiscible droplets in a bifurcating 
channel [64]. Subsequently, the model was generalized through the addition 
of dumbbell-shaped surfactants to model micellization [65] and the behavior 
of ternary amphiphilic mixtures in both two and three dimensions [66,67]. 
Note that since the color current after the collision is always parallel to the 
color gradient, thermal fluctuations of the order parameter are neglected in 
this approach. 

7 

Boundary Conditions and Embedded Objects 
7.1 

Collisional Coupling to Embedded Particles 

A very simple procedure for coupling embedded objects such as colloids or 
polymers to a MPC solvent has been proposed in Ref . [68] . In this approach, 
every colloid particle or monomer in the polymer chain is taken to be a point- 
particle which participates in the SRD collision. If monomer i has mass vrim 
and velocity Wj, the center of mass velocity of the particles in the collision 
cell is 

^ ^ i=i + E^^i w» ^ 

N^m + N^rrim 

where iY,„, is the number of monomers in the collision aAl. A stochastic col- 
lision of the relative velocities of both the solvent particles and embedded 
monomers is then performed in the collision step. This results in an exchange 
of momentum between the solvent and embedded niononua's. The same proce- 
dure can of course be employed for other MPC algorithms, such as MPC- AT. 
The new monomer momenta are then used as initial conditions for a molecular- 
dynamics update of the polymer degrees of freedom during the subsequent 
streaming time step. At. Alternatively, the momentum exchange, Ap, can be 
included as an additional force Ap/At in the molecular-dynamics integration. 
If there are no other interactions between monomers — as might be the case for 
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embedded colloids — these degrees of freedom stream freely during this time 
interval. 

When using this approach, the average mass of solvent particles per cell, 
mNc, should be of the order of the monomer or colloid mass mm (assuming 
one embedded particle per cell) [15]. This corresponds to a neutrally buoyant 
object which responds quickly to the fluid flow but is not kicked around too 
violently. It is also important to note that the average number of monomers 
per cell, (Nm), should be smaller than unity in order to properly resolve 
hydrodynamic interactions between the monomers. On the other hand, the 
average bond length in a semi-flexible polymer or rod-like colloid should also 
not be much larger than the cell size a, in order to capture the anisotropic 
friction of rod-like molecules due to hydrodynamic interactions [69] (which 
leads to a twice as large perpendicular than parallel friction coeSicient for 
long stiff rods [70]), and to avoid an unnecessarily large ratio of the number 
of solvent to solute particles. For a polymer, the average bond length should 
therefore be of the order of a. 

In order to use SRD to model suspended colloids with a radius of order 
1/im in water, this approach would require approximately 60 solvent particles 
per cell in order to match the Peclet number [16]. This is much larger than 
the optimum number (see discussion in Sec. 2.2.1), and the relaxation to the 
Boltzmann distribution is very slow. Because of its simplicity and efficiency, 
this monomer-solvent coupling has been used in many polymer [14,71-74] and 
colloid simulations [15,16,75,76]. 

7.2 

Thermal Boundaries 

In order to accurately resolve the local flow field around a colloid, more ac- 
curate methods have been proposed which exclude fluid-particles from the 
interior of the colloid and mimic slip [19, 77] or no-slip [78] boundary condi- 
tions. The latter procedure is similar to what is known in molecular dynamics 
as a "thermal wall" boundary condition: fluid particles which hit the colloid 
particle are given a new, random velocity drawn from the following probabil- 
ity distributions for the normal velocity component, vn, and the tangential 
component, Vt, 

Pn{vn) = {mvN /ksT) ex.p [—mv%/2kBT) , with vn > , 

Pt{vt) = \frnj2^Kk^ exp {-mv^/2kBT) . (74) 

These probability distributions arc constructed so that the probability distri- 
bution for particles near the wall remains Maxwellian. The probability distri- 
bution, px, for the tangential components of the velocity is Maxwellian, and 
both positive and negative values arc permitted. The normal component must 
be positive, since after scattering at the surface, the particle must move away 
from the wall. The form of pjv is a reflection of the fact that more particles 
with large \vn\ hit the wall per unit time than with small \vn\ [78]. 
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This procedure models a no-slip boundaxy condition at the surface of the 
colloid, and also thermostats the fluid at the boundaries. For many non- 
equilibrium flow conditions, this may not be sufficient, and it may also be 
necessary to thermostat the bulk fluid also (compare Sec. 2.3). It should also 
be noted that Eqs. (74) will be a good approximation only if the radius of 
the embedded objects is much larger than the mean free path A. For smaller 
particles, corrections are needed. 

If a particle hits the surface at time to in the interval between nAt and 
(n + I) At, the correct way to proceed would be to give the particle its new 
velocity and then have it stream the remaining time (n + l)At — to. However, 
such detailed resolution is not necessary. It has been found [16] that good 
results are also obtained using the following simple stochastic procedure. If 
a particle is found to have penetrated the colloid during the streaming step, 
one simply moves it to the boundary and then stream a distance Vnew At e, 
where e is a uniformly distributed random number in the interval [0,1]. 

Another subtlety is worth mentioning. If two colloid particles are very 
close, it can happen that a solvent particle could hit the second colloid after 
scattering off the ffrst, all in the interval At. Naively, one might be tempted to 
simply forbid this from happening or ignore it. However, this would lead to a 
strong depletion- like attractive force between the colloids [16]. This effect can 
be greatly reduced by allowing multiple collisions in which one solvent particle 
is repeatedly scattered off the two colloids. In every collision, momentum is 
transferred to one of the colloids, which pushes the colloids further apart. 
In practice, even allowing for up to ten multiple collisions cannot completely 
cancel the depletion interaction — one needs an additional repulsive force to 
eliminate this unphysical attraction. The same effect can occur when a colloid 
particle is near a wall. 

Careful tests of this thermal coupling have been performed by Padding et 
al. [17,79], who were able to reproduce the correct rotational diffusion of a 
colloid. It should be noted that because the coupling between the solvent par- 
ticles and the surface occurs only through the movement of the fluid particles, 
the coupling is quite weak for small mean free paths. 

7.3 

Coupling using Additional Forces 

Another procedure for coupling an embedded object to the solvent has been 
pursued by Kapral et al. [19,30,80]. They introduce a central repulsive force 
between the solvent particles and the colloid. This force has to be quite strong 
in order to prohibit a large number of solvent particles from penetrating the 
colloid. When implementing this procedure, a small time step 6t is therefore 
required in order to resolve these forces correctly, and a large number of molec- 
ular dynamics time steps are needed during the SRD streaming step. In its 
original form, central forces were used, so that only slip boundary conditions 
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could be modeled. In principle, non-central forces could be used to impose 
no-slip conditions. 

This approach is quite natural and very easy to implement; it docs, how- 
ever, require the use of small time steps and therefore may not be the optimal 
procedure for many applications. 

7.4 

"Ghost" or "Wall" Particles 

One of the first approaches employed to impose a non-slip boundary con- 
dition at an external wall or at a moving object in a MPC solvent was to 
use "ghost" or "wall" particles [36,81]. In other mesoscale methods such as 
Lattice-Boltzmann, no-slip conditions are modeled using the bounce-back rule: 
the velocity of particle is inverted from v to — v when it intersects a wall. For 
planar walls which coincide with the boundaries of the c;ollision c;dls, the same 
procedure can be used in MPC. However, the walls will generally not coincide 
with, or even be parallel to, the cell walls. Furthermore, for small mean free 
paths, where a shift of the cell lattice is required to guarantee Galilean invari- 
ance, partially occupied boundary cells are unavoidable, even in the simplest 
flow geometries. 




Fig. 6. Velocity field of a fluid near a square cylinder in a Poiseuille flow at Reynolds 

number Re — VmaxL/u — 30. The channel width is eight times larger than the 
cylinder size L. A pair of stationary vortices are seen behind the obstacle, as expected 
for Re < 60. Prom Ref. [81]. 

The simple bounce-back rule fails to guarantee no-slip boundary conditions 
in the case of partially filled cells. The following generalization of the bounce- 
back rule has therefore been suggested. For all cells that are cut by walls, 
fill the "wall" part of the cell with a sufficient number of virtual particles in 
order to make the total number of particles equal to M, the average number of 
particles per cell. The velocities of the wall particles are drawn from a Maxwell- 
Boltzmann distribution with zero mean velocity and the same temperature as 
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the fluid. The colUsion step is then carried out using the mean velocity of 
all particles in the cell. Note that since Gaussian random numbers are used, 
and the sum of Gaussian random numbers is also Gaussian-distributed, the 
velocities of the individual wall particles need not to be determined explicitly. 
Instead, the average velocity u can be written as u = +a) /M, where 

a is vector whose components arc Gaussian random numbers with zero mean 
and variance (M — n)kBT. Results for Poiseuille flow obtained using this 
procedure, both with and without cell shifting, were found to be in excellent 
agreement with the correct parabolic profile [36]. Similarly, numerical results 
for the recirculation length, the drag coefficient, and the Strouhal number for 
flows around a circular and square cylinder in two dimensions were shown 
to be in good agreement with experimental results and computational fluid 
dynamics data for a range of Reynolds numbers between Re = 10 and Re = 
130 (see Fig. 6) [36,81]. 

8 

Importance of Angular-Momentum Conservation: Couette 
Flow 

As an example of a situation in which it is important to use an algorithm 
which conserves angular momentum, consider a drop of a highly viscous fluid 
inside a lower-viscosity fluid in circular Couette flow. In order to avoid the 
complications of phase-separating two-component fluids, the high viscosity 
fluid is confined to a radius r < R\hj an impenetrable boundary with re- 
flecting boundary conditions (i.e., the momentum parallel to the boundary is 
conserved in collisions). No-slip boundary conditions between the inner and 
outer fiuids are guaranteed because collision cells reach across the boundary. 
When a torque is applied to the outer circular wall (with no-slip, bounce-back 
boundary conditions) of radius i?2 > a solid-body rotation of both fluids is 
expected. The results of simulations with both MPC-AT-o and MPC-AT+a 
are shown in Fig. 7. While MPC-AT+a reproduces the expected behavior. 
MPC-AT— a produces different angular velocities in the two fluids, with a low 
(high) angular velocity in the fluid of high (low) viscosity [38]. 

The origin of this behavior is that the viscous stress tensor in general has 
symmetric and antisymmetric contributions (see Sec. 4.1.1), 

= X{djV^)dal3 + V {dpVa + d^V^) + f] {dpVa - daV/)) , (75) 

where A is the second viscosity coefficient and f] = pui and fj = pu2 are the 
symmetric and anti-symmetric; components of the viscosity, respectively. The 
last term in Eq. (75) is linear in the vorticity V x v, and does not conserve 
angular momentum. This term therefore vanishes {i.e., fj = 0) when angular 
momentum is c;ons(uvcd. 

The anti-symmetric part of the stress tensor implies an additional torque, 
which becomes relevant when the boundary condition is given by forces. In 
cylindrical coordinates {r,6,z), the azimuthal stress is given by [38] 
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Fig. 7. Azimuthal velocity of binary fluids in a rotating cylinder with J7o = 
0.01(fcsr/moa^)^/^ The viscous fluids with particle mass mi and mo are located 
at r < _Ri and Ri < r < R2, respectively, with _Ri — 5a and R2 = 10a. Symbols 
represent the simulation results of MPC-AT— a with m\/mo — 2 (+) or mi /mo — 5 
(x), and MPC-AT+a for mi /mo = 5 (o). Solid lines represent the analytical results 
for MPC-AT— a at mi /mo = 5. Error bars are smaller than the size of the symbols. 
From Ref. [38]. 



^re = iv + vf-^+2n'^. (76) 
or r 

The first term is the stress of the angular-momentum-conserving fluid, which 
depends on the derivative of the angular velocity f2 = vg/r. The second term 
is an additional stress caused by the lack of angular momentum conservation; 
it is proportional to Q. 

In the case of the phase-separated fluids in circular Couette flow, this 
implies that if both fluids rotate at the same angular velocity, the inner and 
outer stresses do not coincide. Thus, the angular velocity of the inner fluid i7i 
is smaller than the outer one, with vg{r) — ilir for r < Ri and 

ve{r) = Ar + B/r, with A= ^^2 _ ^^2 — - B= _ W^) 

for Ri < r < R2. f^i is then obtained from the stress balance a,t r — Ri, 
i.e., 2fjiQi — (8/3)7/2(^^0 ^ ^1) + 2?72i7i. This calculation reproduces the 
numerical results very well, see Fig. 7. Thus, it is essential to employ an 
+a version of MFC in simulations of multi-phase flows of binary fluids with 
different viscosities. 

There are other situations in which the lack angular momentum conserva- 
tion can cause signiflcant deviations. In Ref. [38], a star polymer with small 
monomer spacing was placed in the middle of a rotating Couette cell. As in the 
previous case, it was observed that the polymer fluid rotated with a smaller 
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angular velocity than the outer fluid. When the angular momentum conser- 
vation was switched on, everything rotated at the same angular velocity, as 
expected. 



9 

MPC without Hydrodynamics 

The importance of hydrodynamic interactions (HI) in complex fluids is gener- 
ally accepted. A standard procedure for determining the influence of HI is to 
investigate the same system with and without HI. In order to compare results, 
however, the two simulations must differ as little as possible — apart from the 
inclusion of HI. A well-known example of this approach is Stokesian dynamics 
simulations (SD), where the original Brownian dynamics (BD) method can 
be extended by including hydrodynamic interactions in the mobility matrix 
by employing the Oseen tensor [12,70]. 

A method for switching off HI in MPC has been proposed in Refs. [24,26]. 
The basic idea is to randomly interchange velocities of all solvent particles 
after each collision step, so that momentum (and energy) axe not conserved 
locally. Hydrodynamic correlations are therefore destroyed, while leaving fric- 
tion coefficients and ffuid self-diffusion coefficients largely unaffected. Since 
this approach requires the same numerical effort as the original MPC algo- 
rithm, a more efficient method has been suggested recently in Ref. [25] . If the 
velocities of the solvent particles are not correlated, it is no longer necessary 
to follow their trajectories. In a random solvent, the solvent-solute interac- 
tion in the collision step can thus be replaced by the interaction with a heat 
bath. This strategy is related to that proposed in Ref. [36] to model no-slip 
boundary conditions of solvent particles at a planar wall, compare Sec. 7.4. 
Since the positions of the solvent particles within a cell are not required in 
the collision step, no explicit particles have to be considered. Instead, each 
monomer is coupled with an effective solvent momentum P which is directly 
chosen from a Maxwell-Boltzmann distribution of variance mMkBT and a 
mean given by the average momentum of the fluid field which is zero at rest, 
or {mMjr'y, 0,0) in the case of an imposed shear flow. The total center-of-mass 
velocity, which is used in the collision step, is then given by [25] 

rrimyi + P 

m,M + m,m 

where is the mass of the solute particle. The solute trajectory is then 
determined using MD, and the interaction with the solvent is performed every 

collision time At. 

The random MPC solvent therefore has similar properties to the MPC 
solvent, except that there are no HI. The relevant parameters in both methods 

are the average number of particles per cell, M, the rotation angle a, and the 
collision time At which can be chosen to be the same. For small values of the 
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density (M < 5), fluctuation effects have been noticed [26] and could also be 
included in the random MPC solvent by a Poisson-distributed density. The 
velocity autocorrelation functions [15] of a random MPC solvent show a simple 
exponentially decay, which implies some differences in the solvent diffusion 
coefficients. Other transport coefficients such as the viscosity depend on HI 
only weakly [57] and consequently are expected to be essentially identical in 
both solvents. 

10 

Applications to Colloid and Polymer Dynamics 

The relevance of hydrodynamic interactions for the dynamics of complex 
fluids — such as dilute or semidilute polymer solutions, colloid suspensions, 
and microcmulsions — is well known [12, 70]. From the simulation point of 
view, however, these systems are difficult to study because of the large gap in 
length- and time-scales between solute and solvent dynamics. One possibility 
for investigating complex fluids is the straightforward application of molec- 
ular dynamics simulations (MD), in which the fluid is course-grained and 
represented by Lennard- Jones particles. Such simulations provide valuable in- 
sight into polymer dynamics [83-87]. Similarly, mesoscale algorithms such as 
Lattice-Boltzmann and dissipative particle dynamics have been widely used 
for modeling of colloidal and polymeric systems [88-92]. 

Solute molecules, e.g., polymers, arc typically composed of a large number 
of individual particles, whose interactions are described by a force-fleld. As 
discussed in Sec. 7, the particle-based character of the MPC solvent allows 
for an easy and controlled coupling between the solvent and solute particles. 
Hybrid simulations combining MPC and molecular dynamics simulations are 
therefore easy to implement. Results of such hybrid simulations are discussed 
in the following. 

10.1 
Colloids 

Many applications in chemical engineering, geology, and biology involve sys- 
tems of particles immersed in a liquid or gas flow. Examples include sedimen- 
tation processes, liquid-solid fluidized beds, and flocculation in suspensions. 

Long-ranged solvent-mediated hydrodynamic interactions have a profound 
effect on the non-equilibrium properties of colloidal suspensions, and the 
many-body hydrodynamic backflow effect makes it difficult to answer even 
relatively simple questions such what happens when a collection of parti- 
cles sediments through a viscous fluid. Batchelor [93] calculated the lowest- 
order volume fraction correction to the average sedimentation velocity, Vg = 
v'^{l — 6.55 (j!)), of hard spheres of hydrodynamic radius Rh where is the 
sedimentation velocity of a single sphere. Due to the complicated interplay 
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between short-ranged contact forces and fong-ranged HI, it is hard to extend 
this result to the high volume fraction suspensions of interest for ceramics 
and soil mechanics. An additional complication is that the Brownian motion 
of solute particles in water cannot be neglected if they are smaller than 1/um 
in diameter. 

The dimcnsionlcss Pcclct number Pe = v'^Rh/D, where D is the self- 
diffusion coefficient of the suspended particles, measures the relative strength 
of HI and thermal motion. Most studies of sedimentation have focused on 
the limit of infinite Pcclct number, where Brownian forces arc negligible. For 
example, Ladd [94] employed a Lattice-Boltzmann method (LB), and Hoefler 
and Schwarzer [95] used a marker-and-cell Navier-Stokes solver to simulate 
such non-Brownian suspensions. The main difficulty with such algorithms is 
the solid-fluid coupling which can be very tricky: in LB simulations, special 
"boundary nodes" were inserted on the colloid surface, while in Ref. [95], the 
coupling was mediated by inertia-less markers which are connected to the 
colloid by stiff springs and swim in the fluid, effectively dragging the colloid, 
but also exerting a force on the fluid. Several methods for coupling embedded 
particles to an MFC solvent were discussed in Sec. 7. 

Using the force-based solvent-colloid coupling described in Sec. 7.3, Fadding 
and Louis [96] investigated the importance of HI during sedimentation at small 
Fcclct numbers. Surprisingly, they found that the sedimentation velocity docs 
not change if the Feclet number is varied between 0. 1 and 15 for a range of vol- 
ume fractions. For small volume fractions, the numerical results agree with the 
Batchclor law; for intermediate (j) they are consistent with the semi-empirical 
Richardson-Zaki law, Vg = Vg{l — </))", n = 6.55. Even better agreement was 
found with theoretical predictions by Hayakawa and Ichiki [97, 98] , who took 
higher order HI into account. Furcly hydrodynamic arguments arc therefore 
still valid in an average sense at low Pe, i.e., for strong Brownian motion and 
relatively weak HI. This also means that pure Brownian simulations without 
HI. which lead to Vg = (1 — (^), strongly underestimate the effect of backfiow. 

On the other hand, it is known that the velocity autocorrelation func- 
tion of a colloidal particle embedded in a fluctuating liquid at equilibrium 
exhibits a hydrodynamic long-time tail, (u(t)f(O)) ^ t~'^/'^, where d is the 
spatial dimension [99]. These tails have been measured earlier for point-like 
SRD particles in two [21,27] and three [15] spatial dimensions, and found to 
be in quantitative agreement with analytic predictions, with no adjustable 
parameters. It is therefore not surprising that good agreement was also ob- 
tained for embedded colloids [96]. MFC therefore correctly describes two of 
the most important effects in colloidal suspensions, thermal fluctuations and 
hydrodynamic interactions. 

In a series of papers, Hecht et al. [16, 76, 100] used hybrid SRD-molecular 
dynamics simulations to investigate a technologically important colloidal 
system — yl?203-particles of diameter 0.5/im (which is often used in ceram- 
ics) suspended in water — ^with additional colloid-colloid interactions. These 
colloids usually carry a charge which, by forming an electric double layer with 
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ions in water, results in a screened electrostatic repulsion. The interaction can 
be approximated by the Derjaguin-Landau-Verwey-Overbeek (DLVO) the- 
ory [101, 102]. The resulting potential contains a repulsive Dcbyc-Huckcl con- 
tributions, Vel ~ exp(— K;[r — cJ])/r, where d is the particle diameter, k is 
the inverse screening length, and r is the distance of the particle centers. The 
second part of the DLVO-potential is a short-range van der Waals attraction. 



which turns out to be important at the high volume fractions (</> > 20%) and 
high salt concentrations of interest. Ah is the Hamaker constant which in- 
volves the polarizability of the particles and the solvent. DLVO theory makes 
the assumption of linear polarizability and is valid only at larger distances. It 
therefore does not include the so-called primary potential minimum at particle 
contact, which is observed experimentally and is about SOkgT deep. Because 
of this potential minimum, colloids which come in contact rarely become free 
again. In order to ensure numerical stability for reasonable values of the time 
step, this minimum was modeled by an additional parabolic potential with 
depth of order QUbT. The particle Reynolds number of the real system is 
very small, of order 10~^ to 10"''. Since it would be too time-consuming to 
model this Reynolds number, the simulations were performed at Re w 0.02, 
which still ensures that the contribution of momentum convection is negligi- 
ble compared to that of momentum diffusion. However, due to the remaining 
incrtial effects and the non-zero time step, it was still possible that particles 
partially overlapped in the simulation. This overlap was penalized by an ad- 
ditional potential, frequently used in simulations of granular matter, given by 
a Hertz-law, 



SRD correctly describes long-range HI, but it can only resolve hydrody- 
namic interactions on scales larger than both the mean free path A and the cell 
size a. In a typical simulation with about 1000 colloid particles, a relatively 
small colloid diameter of about four lattice units was chosen for computational 
efficiency. This means that HI are not fully resolved at interpaxticle distances 
comparable to the colloid diameter, and lubrication forces have to be inserted 
by hand. Only the most divergent mode, the so-called squeezing mode, was 
used. Flub ~ Vrei/rreh where r^ei and Vrei are the relative distance and ve- 
locity of two colloids, respectively. This system of interacting ^Z203-particles 
was simulated in order to study the dependence of the suspension's viscosity 
and structure on shear rate, pH, ionic strength and volume fraction. The re- 
sulting stability diagram of the suspension as a function of ionic strength and 
pH value is shown in Fig. 8 (plotted at zero shear) [76]. The pH controls the 
surface charge density which, in turn, affects the electrostatic interactions be- 
tween the colloids. Increasing the ionic strength, experimentally achieved by 
"adding salt" , decreases the screening length 1 /k, so that the attractive forces 
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become more important; the particles start forming clusters. Three different 
states are observed: (i) a clustered regime, where particles aggregate when van 
dcr Waals attractions dominate, (ii) a suspended regime where particles are 
distributed homogeneously and can move freely — corresponding to a stable 
suspension favored when electrostatic repulsion prevents clustering but is not 
strong enough to induce order. At very strong Coulomb repulsion the repul- 
sive regime (iii) occurs, where the mobility of the particles is restricted, and 
particles arrange in local order which maximizes nearest neighbor distances. 

The location of the phase boundaries in Fig. 8 depends on the shear rate. In 
the clustered phase, shear leads to a breakup of clusters, and for the shear rate 
7 = 1000 s~^, there are many small clusters which behave like single particles. 
In the regime where the particles are slightly clustered, or suspended, shear 
thinning is observed. Shear thinning is more pronounced in the slightly clus- 
tered state, because shear tends to reduce cluster size. Reasonable agreement 
with experiments was achieved, and discrepancies were attributed to poly- 
dispersity and the manner in which lubrication forces were approximated, as 
well as uncertainties how the pH and ionic strength enter the model force 
parameters. 




Fig. 8. Phase diagram of a colloidal suspension (plotted at zero shear and volume 
fraction 4> = 35%) in the ionic-strength-pH plane depicting three regions: a clustered 
region, a suspended regime, and a repulsive structure. Prom Ref. [76]. 
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In the simulations of Hecht et al. [16], the simple collisional coupling proce- 
dure described in Sec. 7.1 was used. This means that the colloids were treated 
as point particles, and solvent particles could flow right through them. Hydro- 
dynamic interactions were therefore only resolved in an average sense, which 
is acceptable for studies of the general properties of an ensemble of many 
colloids. The heat from viscous heating was removed using the stochastic 
thermostat described in Sec. 2.3. 

Various methods for modeling no-slip boundary conditions at colloid 
surfaces — such as the thermal wall coupling described in Sec. 7.2 were sys- 
tematically investigated in Ref. [79]. No-slip boundary conditions are impor- 
tant, since colloids are typically not completely spherical or smooth, and the 
solvent molecules also transfer angular momentum to the colloid. Using the 
SRD algorithm without angular momentum conservation, it was found that 
the rotational friction coefficient was larger than predicted by Enskog-theory 
when the ghost-particle coupling was used [82]. On the other hand, in a de- 
tailed study of the translational and rotational velocity autocorrelation func- 
tion of a sphere coupled to the solvent by the thermal-wall boundary condition, 
quantitative agreement with Enskog theory was observed at short times, and 
with mode-coupling theory at long times. However, it was also noticed that 
for small particles, the Enskog and hydrodynamic contributions to the friction 
coefficients were not clearly separated. Specifically, mapping the system to a 
density matched colloid in water, it appeared that the Enskog and the hydro- 
dynamic contributions are equal at a particle radius of 6 nm for translation 
and 35.4 nm for rotation; even for a particle radius of 100 nm, the Enskog 
contribution to the friction is still of order 30% and cannot be ignored. 

In order to clarify the detailed character of the hydrodynamic interactions 
between colloids in SRD, Lee and Kapral [103] numerically evaluated the fixed- 
particle friction tensor for two nano-spheres embedded in an SRD solvent. 
They found that for intercolloidal spacings less than 1.2 d, where d is the 
colloid diameter, the measured friction coefficients start to deviate from the 
expected theoretical curve. The reader is referred to the review by Kapral [30] 
for more details. 

10.2 

Polymer Dynamics 

The dynamical behavior of macromolecules in solution is strongly affected or 
even dominated by hydrodynamic interactions [70,104,105]. From a theoret- 
ical point of view, scaling relations predicted by the Zimm model for, e.g., 
the dependencies of dynamical quantities on the length of the polymer are, in 
general, accepted and confirmed [106]. Recent advances in experimental single- 
molecule techniques provide insight into the dynamics of individual polymers, 
and raise the need for a quantitative theoretical description in order to deter- 
mine molecular parameters such as diffusion coefficients and relaxation times. 
Mesoscale hydrodynamic simulations can be used to verify the validity of the- 
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oretical models. Even more, such simulations are especially valuable when 
analytical methods fail, as for more complicated molecules such as polymer 
brushes, stars, ultrasoft colloids, or scmidihitc solutions, where hydrodynamic 
interactions are screened to a certain degree. Here, mesoscale simulations still 
provide a full characterization of the polymer dynamics. 

Wo will focus on the dynamics of polymer chains in dilute solution. In order 
to compare simulation results with theory — in particular the Zimm approach 
[70, 107] — and scaling predictions, we address the dynamics of Gaussian as 
well as self-avoiding polymers. 



10.2.1 

Simulation Method and Model 

Polymer molecules are composed of a large number of equal repeat units called 
monomers. To account for the generic features of polymers, such as their con- 
formational freedom, no detailed modeling of the basic units is necessary. A 
coarse-grained description often suffices, where several monomers are com- 
prised in an effective particle. Adopting such an approach, a polymer c;haiii 
is introduced into the MPC solvent by adding point particles, each of 
mass rrim, which are connected linearly by bonds. Two different models are 
considered, a Gaussian polymer and a polymer with cxcluded-volume (EV) 
interactions. Correspondingly, the following potentials are applied: 

(i) Gaussian chain: The monomers, with the positions r, {i = 1, . . . , Nm), are 
connected by the harmonic potential 

t^« = ^ E'('-«+i-rO% (81) 

i=l 

with zero mean bond length, and b the root-mean-square bond length. Here, 
the various monomers freely penetrated each other. This simplification allows 
for an analjdiical treatment of the chain dynamics as in the Zimm model 

[70,107]. 

(ii) Excluded-volume chain: The monomers are connected by the harmonic 
potential 

f^B = f E {\ri+i-ri\-bf, (82) 

i=l 

with mean bond length b. The force constant k is chosen such that the fluc- 
tuations of the bond lengths are on the order of a percent of the mean bond 
length. In addition, non-bonded monomers interact via the repulsive, trun- 
cated Lennard-Jones potential 



Ulj 



4e 



, 12 



0, otherwise 



(83) 
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The excluded volume leads to swelling of the polymer structure compared to a 
Gaussian chain, which is difficult to fully account for in analytical calculations 
[73]. 

The dynamics of the chain monomers is determined by Newtons' equations 
of motion between the collisions with the solvent. These equations are inte- 
grated using the velocity Vcrlct algorithm with the time stop Atp. The latter 
is typically smaller than the collision time At. The monomer-solvent interac- 
tion is taken into account by inclusion of the monomer of mass = pm 
in the collision step [68,73], compare Sec. 7.1. Alternatively, a Lcnnard- Jones 
potential can be used to account for the monomer-MPC particle interaction, 
where a MPC particle is of zero interaction range [19,108]. 

We scale length and time according to x = x/a and i = tyjh/^rjma? , 
which corresponds to the choice ksT = 1, m = 1, and a = 1. The mean free 
path of a fluid particle At^/ksT /m is then given by A = At. In addition, we 
set b = a, (7 = a, and e/kBT = 1. 

The equilibrium properties of a polymer are not affected by hydrodynamic 
interactions. Indeed, the results for various equilibrium quantities — such as 
the radius of gyration of MPC simulation arc in excellent agreement with 
the results of molecular dynamics of Monte Carlo simulations without explicit 
solvent [73]. 

Simulations of Gaussian chains, i.e., polymers with the bond potential 
(81), can be compared with analytical calculations based on the Zimm ap- 
proach [70,107]. Note, however, that the simulations are not performed in the 
Zimm model. The Zimm approach relics on the prcavcraging approximation 
of hydrodynamic interactions, whereas the simulations talce into account the 
configurational dependence of the hydrodynamic interactions, and therefore 
hydrodynamic fluctuations. Hence, the comparison can serve as a test of the 
validity of the approximations employed in the Zimm approach. 

The Zimm model rests upon the Langevin equation for over-damped mo- 
tion of the monomers, i.e., it applies for times larger than the Brownian time 
scale tb ^ nim/C, where ( is Stokes' friction coefficient [12]. On such time 
scales, velocity correlation functions have decayed to zero and the monomer 
momenta arc in equilibrium with the solvent. Moreover, hydrodynamic in- 
teractions between the various parts of the polymer are assumed to propa- 
gate instantaneously. This is not the case in our simulations. First of all, the 
monomer inertia term is taken into account, which implies non-zero velocity 
autocorrelation functions. Secondly, the hydrodynamic interactions build up 
gradually. The center-of-mass velocity autocorrelation function displayed in 
Fig. 9 reflects these aspects. The correlation function exhibits a long-time tail, 
which decays as {v cm{t)'v cm{0)) ~ t~^/^ on larger time scales. The algebraic 
decay is associated with a coupling between the motion of the polymer and 
the hydrodynamic modes of the fluid [99, 109, 110]. A scaling of time with 
the diffusion coefficient D shows that the correlation function is a universal 
function of Dt. This is in agreement with results of DPD simulations of dilute 
polymer systems [92]. 
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Fig. 9. Center-of-mass velocity autocorrelation functions for Gaussian polymers of 
length Nm = 20, Nm = 40, and A = 0.1 as a function of Dt. The solid line is 
proportional to {Dt)~^^^. From Ref. [73]. 



The polymer center-of-mass diffusion coefficient follows either via the 
Green-Kubo relation from the velocity autocorrelation function or by the Ein- 
stein relation from the center-of-mass mean square displacement. According 
to the Kirkwood formula [104, 105, 111] 



jj(K) ^ ^ 



1 



where the hydrodynamic radius Rh is defined as 




(84) 



(85) 



and the prime indicates that the term with j = i has to be left out in the 
summation. The diffusion coefficient is composed of the local friction contri- 
bution Dq/Nj^, where Dq is the diffusion coefficient of a single monomer in 
the same solvent, and the hydrodynamic contribution. 

Simulation results for the hydrodynamic contribution, Dh = D — Do/Nm, 
to the diffusion coefficient are plotted in Fig. 10 as a function of the hydrody- 
namic radius (85). In the limit Nm ^ 1, the diffusion coefficient D is domi- 

— 1/2 

nated by the hydrodynamic contribution Dh, since Dh ~ Nm ■ For shorter 
chains, D^/Nm cannot be neglected, and therefore has to be subtracted in 
order to extract the scaling behavior of Dh- The hydrodynamic part of the 
diffusion coefficient Dh exhibits the dependence predicted by the Kirkwood 
formula and the Zimm theory, i.e., Dh ~ ^/Rh- The finite-size corrections to 
D show a dependence D = D^o — const. /L on the size L of a periodic system. 
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Fig. 10. Dependence of the hydrodynamic part of the diffusion coefficient, Dh = 
D — Do/Nm, on the hydrodynamic radius for Gaussian chains of lengths Nm = 5, 
10, 20, 40, 80, and 160 (right to left). The mean free path is A = 0.1. Prom Ref. [73]. 

in agreement with previous studies [68,91,112]. Simulations for various sys- 
tem sizes for polymers of lengths A^^, = 10, 20, and 40 allow an extrapolation 
to infinite system size, which yields Dq/ y^ksTa"^ /m « 1.7 x 10~^, in good 
agreement with the diffusion coefficient of a monomer in the same solvent. 
The values of £>oo are about 30% larger than the finite-system-size values pre- 
sented in Fig. 10. Similarly the diffusion coefficient for a polymer chain with 
exchidcd volume interactions displays the dependence Dh ~ ^/Rh [73]. 

The Kirkwood formula neglects hydrodynamic fluctuations and is thus 
identical with the preaveraging result of the Zimm approach. When only the 
hydrodynamic part is considered, the Zimm model yields the diffusion coeffi- 
cient 



MFC simulations for polymers of length Nm = 40 yield Dz / y^ksTa'^ /m = 
0.003. This value agrees with the numerical value for an infinite system, 
Dn/^/kBTa^/m = 0.0027, within 10%. The MFC simulations yield a dif- 
fusion coefficient smaller than D^^\ in agreement with previous studies pre- 
sented in Refs. [70,111, 113]. Note that the experimental values are also smaller 
by about 15% than those predicted by the Zimm approach [70, 114, 115]. 

To further characterize the internal dynamics of the molecular chain, a 
mode analysis in terms of the eigenfunctions of the discrete Rouse model [70, 
116] has been performed. The mode amplitudes Xp are calculated according 
to 



Dz = 0.192 



(86) 
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Fig. 11. Correlation functions of the Rouse-mode amplitudes for the modes p — l — i 
of Gaussian polymers. The chain lengths are Nm = 20 (right) and Nm = 40 (left). 
From Ref. [73]. 



Due to hydrodynamic interactions, Rouse modes are no longer eigenfunctions 
of the chain molecule. However, within the Zimm theory, they are reasonable 
approximations and the autocorrelation functions of the mode amplitudes 
decay exponentially, i.e., 

(Xp(i)Xp(O)) = (xl) exp i-t/Tp) . (88) 

For the Rouse model, the relaxation times Tp depend on chain length and mode 
number according to Tp ~ 1/sin'^ (p7r/2iVm), whereas for the Zimm model the 
dependence 

Tp ^ (p/iV™)^/V sin' ipn/2Nm) (89) 

is obtained. The extra contribution \/p/N.m follows from the eigenfunction 
representation of the preaveraged hydrodynamic tensor, under the assumption 
that its off-diagonal elements do not significantly contribute to the relaxation 
behavior. 

In Fig. 11, the autocorrelation functions for the mode amplitudes are 
shown for the mean free path A — 0.1. Within the accuracy of the simula- 
tions, the correlation functions decay exponentially and exhibit the scaling 
behavior predicted by the Zimm model. Hence, for the small mean free path, 
hydrodynamic interactions are taken into account correctly. This is no longer 
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Fig. 12. Dependence of the longest relaxation time n on the radius of gyration for 
Gaussian chains of the lengths given in Fig. 10. Prom Ref. [73] . 

true for the large mean free path, A = 2. In this scaHng behavior be- 

tween that predicted by the Rouse and Zimm models is observed. This implies 
that hydrodynamic interactions are present, but are not fully developed or are 
small compared to the local friction of the monomers. We obtain pure Rouse 
behavior for a system without solvent by simply rotating the velocities of the 
individual monomers [73]. 

The dependence of the longest relaxation time on the radius of gyration 
is displayed in Fig. 12 for A = 0.1. The scaling behavior ti ^ Rg is in very 
good agreement with the predictions of the Zimm theory. We even find al- 
most quantitative agreement; the relaxation time of the p = 1 mode of our 
simulations is approximately 30 % larger than the Zimm value [70]. 

The scaling behavior of equilibrium properties of single polymers with 
excluded- volume interactions has been studied extensively [70, 117-120]. It 
has been found that even very short chains already follow the scaling be- 
havior expected for much longer chains. In particular, the radius of gyration 
increases like Rg ~ ./V^ with the number of monomers, and the static struc- 
ture factor S{q) exhibits a scaling regime for 2tt /Rg <C q ^ 27r/(T, with a 
g-i/'' decay as a function of the scattering vector q and the exponent v w 0.6. 
For the interaction potentials (82), (83) with the parameters h = a = a, 
e/ksT = 1, the exponent v « 0.62 is obtained from the chain-length depen- 
dence of the radius of gyration, the mean square end-to-end distance, as well 
as the dependence of the static structure factor [73]. 

An analysis of the intramolecular dynamics in terms of the Rouse modes 
yields non-exponentially decaying autocorrelation functions of the mode am- 
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Fig. 13. Correlation functions of the Rouse-mode amplitudes for various modes as a 
function of the scaled time tp°' for polymers with excluded volume interactions. The 
chain lengths are Nm = 20 (left) and Nm = 40 (right). The calculated correlations 
where fitted by Ap exp(— t/r^) and have been divided by Ap. The scaling exponents 
of the mode numbers are a = 1.93 (A'^^, = 20) and a = 1.85 (A'^^ = 40), respectively. 
From Ref. [73]. 



plitudes. At very short times, a fast decay is found, which turns into a slower 
exponential decay which is well fitted by Ap exp{—t /Tp) , see Fig. 13. Within 
the accuracy of these calculations, the correlation functions exhibit universal 
behavior. Zimm theory predicts the dependence Tp p~^^ for the relax- 
ation times on the mode number for polymers with excluded-volume interac- 
tions [70]. With V = 0.62, the exponent a for the polymer of length N„i = 40 
is found to be in excellent agreement with the theoretical prediction. The 
exponent for the polymers with — 20 is slightly larger. 

Zimm theory predicts that the dynamic structure factor, which is defined 

by 

'^('l' = ]^ E E (e^P - ^^(0)])) , (90) 

™ i=l j = l 

scales as [70] 

S{q,t) = S{q,0)f{q^t) (91) 

with a = 3 for qRo ^ 1, independent of the solvent conditions (0 or good sol- 
vent). To extract the scaling relation for the intramolecular dynamics, which 
corresponds to the prediction (91), we resort to the following considerations. 
As is well known, the dynamic structure factor for a Gaussian distribution of 
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Fig. 14. Normalized dynamic structure factor ^(q, t)/(5'(q, 0) exp (—Dq^f)) of poly- 
mers with excluded volume interactions for Nm = 40 and various g-values in the 
range 0.7 < ga < 2 as a function of g^t^/^. From Ref. [73]. 



the differences ri{t)—rj{0) and a linear equation of motion is given by [70, 121] 

S{ci,t) = 5(q,0)e-^^^*— 5]^exp(V((r:W -r;-(0))^)/6) , (92) 

™ 1=1 j=i 

where Dq^t accounts for the center-of-mass dynamics and denotes the 
position of monomer i in the center-of-mass reference frame. Therefore, 
in order to obtain the dynamics in the center-of-mass reference frame, we 
plot 5(q, t)/{S{q, 0) exp {—Dq^tj). The simulation results for the polymer of 
length iV,„ = 40, shown in Fig. 14, confirm the predicted scaling behavior. 
Thus, MPC-MD hybrid simulations are very well suited to study the dynam- 
ics of even short polymers in dilute solution. 

As mentioned above, the structure of a polymer depends on the nature of 
the solvent. In good solvent excluded volume interactions lead to expanded 
conformations and under bad solvent conditions the polymer forms a dense 
coil. In a number of simulations the influence of hydrodynamic interactions on 
the transition from an extended to a collapsed state has been studied, when 
the solvent quality is abruptly changed. Both, molecular dynamics simula- 
tions with an explicit solvent [122] as well as MFC simulations [108, 123, 124] 
yield a significant different dynamics in the presence of hydrodynamic inter- 
actions. Specifically, the collapse is faster, with a much weaker dependence of 
the characteristic time on polymer length [108, 124], and the folding path is 
altered. 

Similarly, a strong infiuence of hydrodynamic interactions has been found 
on the polymer translocation dynamics through a small hole in a wall [125] or 
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in polymer packing in a virus capsid [126, 127]. Cooperative backflow effects 
lead to a rather sharp distribution of translocation times with a peak at rela- 
tively short times. The fluid flow field, which is created as a monomer moves 
through the hole, guides following monomers to move in the same direction. 

10.3 

Polymers in Flow Fields 

Simulations of a MPC fluid confined between surfaces and exposed to a con- 
stant external force yield the expected parabohc velocity profile for appropri- 
ate boundary conditions [31,36, 128]. The abihty of MPC to account for the 
flow behavior of mesoscale objects, such as polymers, under non-equilibrium 
conditions has been demonstrate for a number of systems. Rod-like colloids 
in shear flow exhibit flow induced alignment [72]. The various diagonal com- 
ponents of the radius of gyration tensor exhibit a qualitative and quantitative 
different behavior. Due to the orientation, the component in the flow direc- 
tion increases with increasing Peclet number larger than united and saturates 
at large shear rates because of finite size effects. The transverse components 
decrease with shear rate, where the component in the gradient direction is 
reduced to a greater extent. The rod rotational velocity in the shear plane 
shows two distinct regimes. For Peclet numbers much smaller than unity, the 
rotational velocity increases linearly with the shear rate, because the system 
is isotropic. At Peclet numbers much larger than unity, the shear induced 
anisotropies lead to a slower increase of the rotational velocity with the shear 
rate [72]. 

The simulations of a tethered polymer in a Poiseuille flow [74] yield a series 
of morphological transitions from sphere to deformed sphere to trumpet to 
stem and flower to rod, similar to theoretically predicted structures [129- 
131]. The crossovers between the various regimes occur at flow rates close to 
the theoretical estimates for a similar system. Moreover, the simulations in 
Ref. [74] show that backflow effects lead to an effective increase in viscosity, 
which is attributed to the fluctuations of the free polymer end rather than its 
shape. 

The conformational, structural, and transport properties of free flexible 
polymers in microchannel flow have been studied in Refs. [128, 132] by hybrid 
MPC-molecular dynamics simulations. These simulations confirm the cross 
streamline migration of the molecules as previously observed in Refs. [133- 
138]. In addition, various other polymer properties are addressed in Ref. [132]. 

All these hybrid simulations confirm that MPC is an excellent method to 
study the non-equilibrium behavior of polymers in flow flelds. In the next sec- 
tion, we will provide a more detailed example for a more complicated object, 
namely an ultrasoft colloid in shear flow. 



56 



G. Gompper, T. Ihle, D.M. KroU, and R.G. Winkler 




Fig. 15. Fluid flow lines in the flow-gradient plane of the star polymer's center of 
mass reference frame for / — 10 (top) and f — 50 (bottom) arms, both with L/ = 20 
monomers per arm and an applied shear field with Wi = -yr = 22. From Ref. [25]. 



10.4 

Ultra-soft Colloids in Shear Flow 

Star polymers present a special macromolecular architecture, in which sev- 
eral linear polymers of identical length are linked together by one of their 
ends at a common center. This structure is particularly interesting because it 
allows for an almost continuous change of properties from that of a flexible 
linear polymer to a spherical colloidal particle with very soft interactions. Star 
polymers are therefore also often called ultrasoft colloids. The properties of 
star polymers in and close to equilibrium have been studied intensively, both 
theoretically [139, 140] and experimentally [141]. A star polymer is a ultrasoft 
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Fig. 16. Normalized component Gxx of the average gyration tensor as a function 
of the Weissenberg number Wi, for star polymers of / = 50 arms and the arm 
lengths Lf = 10, 20, and 40 monomers. Power-law behaviors with quadratic and 
linear dependencies on Wi are indicated by lines. The inset shows all three diagonal 
components of the gyration tensor (for I// = 20). From Ref. [142]. 



colloid, where the core extension is very small compared to the length of an 
arm. By anchoring polymers on the surface of a hard colloid, the softness 
can continuously be changed from ultrasoft to hard by increasing the ratio 
between the core and shell radius at the expense of the thickness of the soft 
polymer corona. Moreover, star polymers have certain features in common 
with vesicles and droplets. Although their shell can be softer than that of the 
other objects, the dense packing of the monomers will lead to a cooperative 
dynamical behavior resembling that of vesicles or droplets [142]. 

Vesicles and droplets encompass fluid which is not exchanged with the 
surrounding. In contrast, for star-like molecules fluid is free to penetrate into 
the molecule and internal fluid is exchanged with the surrounding in the course 
of time. This intimate coupling of the star-polymer dynamics and the fluid 
flow leads to a strong modification of the flow behavior at and next to the 
ultrasoft colloid particular in non-equilibrium systems. 

In the following, we will discuss a few aspects in the behavior of star 
polymers in shear flow as a function of arm number /, arm length Lf, and 
shear rate 7. The polymer model is the same as described in Sec. 10.2.1, 
where the chain connectivity is determined by the bond potential (82) and 
the excluded-volume interaction is described by the Lennard- Jones potential 
(83). A star polymer of functionality / is modeled as / linear polymer chains 
of Lf monomers each, with one of their ends linked to a central particle. 
Linear polymer molecules are a special case of star polymers with function- 
ality / = 2. Lees-Edwards boundary conditions [42] are employed in order to 
impose a linear velocity profile {vx,Vy, Vz) = {i^y, 0, 0) in the fluid in the ab- 
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sence of a polymer. For small shear rates, the conformations of star polymers 
remain essentially unchanged compared to the equilibrium state. Only when 
the shear rate exceeds a characteristic value, a structural anisotropy as well 
as an alignment is induced by the flow (cf. Fig. 15). The shear rate depen- 
dent quantities are typically presented in terms of the Weissenberg number 
Wi = 7T rather than the shear rate itself, where t is the longest characteristic 
relaxation time of the considered system. For the star polymers, the best data 
collapse for stars of various arm lengths is found when the relaxation time 
r = rjb^L'j/kBT is used [142]. Remarkably, there is essentially no dependence 
on the functionality. Within the range of investigated star sizes, this relax- 
ation time has to be considered as consistent with the prediction for the blob 
model of Ref. [143], where r ^ Lj.-^/"-^ for the Flory exponent z/ = 0.6. 

In Fig. 15, typical star conformations are shown which indicate the align- 
ment and induced anisotropy in the flow. Moreover, the figure reveals the 
intimate coupling of the polymer dynamics and the emerging fluid flow field. 
In the region, where the fiuid coexists with the star polymer, the externally im- 
posed fiow field is strongly screened and the fiuid velocity is no longer aligned 
with the shear flow direction, but rotates around the polymer center of mass. 
The fluid stream lines are calculated by integration of the coarse-grained fluid 
velocity field. Outside the region covered by the star polymer, the fiuid adapts 
to the central rotation by generating two counter-rotating vortices, and cor- 
respondingly two stagnation points of vanishing fiuid velocity [25]. 

The fiuid fiow in the vicinity of the star polymer is distinctively different 
from that of a sphere but resembles the fiow around an ellipsoid [144]. In 
contrast to the latter, the fiuid penetrates into the area covered by the star 
polymer. While the fiuid in the core of the star rotates together with the 
polymer, the fluid in the corona follows the external flow to a certain extent. 

A convenient quantity to characterize the structural properties and align- 
ment of polymers in fiow is the average gyration tensor, which is defined as 

G„M7) = i^E<<X/3) . (93) 

"* i=l 

where Nm = fLf + 1 is the total number of monomers, r- ^ is the position of 

monomer i relative to the polymer center of mass, and a G {x, y, z}. The av- 
erage gyration tensor is directly accessible in scattering experiments. Its diag- 
onal components Gaa(7) are the squared radii of gyration of the star polymer 
along the axes of the reference frame. In the absence of flow, scaling consid- 
erations predict [139] G^^{0) = Gyy{0) = G^^{0) = i?2(0)/3 ~ Lf f-" . 

The diagonal components Gaa of the average gyration tensor are shown 
in Fig. 16 as a func;tion of the Weissenberg number for various functionalities 
and axm lengths. We find that the extension of a stax increases with increasing 
shear rate in the shear direction {x), decreases in the gradient direction (y), 
and is almost independent of Wi in the vorticity direction (z). The deviation 
from spherical symmetry exhibits a Wi^ power-law dependence for small shear 
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Fig. 17. Orientational resistance mo as a function of the Weissenberg number Wi 
for star polymers with functionahties / = 2 (circles), / = 15 (triangles), and / = 50 
(squares), and different arm lengths indicated in the figure. Lines correspond to the 
power law mc ~ mG(/)Wi'''^^. The inset shows that the amplitude also follows a 
power-law behavior with rha{f) ~ Z"'^^- From Ref. [142]. 

rates for all functionalities. A similar behavior has been found for rod- like 
colloids [72] (due to the increasing alignment with the flow direction) and for 
linear polymers [70]. However, for stars of not too small functionality, a new 
scaling regime appears, where the deformation seems to scales linearly with 
the Weissenberg number. For large shear rates, finite-size effects appear due 
to the finite monomer number. These effects emerge when the arms are nearly 
stretched, and therefore occur at higher Weissenberg numbers for larger arm 
lengths. 

The average ffow alignment of a (star) polymer can be characterized by 
the orientation angle Xd which is the angle between the eigenvector of the 
gyration tensor with the largest eigenvalue and the flow direction. It follows 
straightforwardly [87] from the simulation data via 



where the right-hand-side of the equation deflnes the orientation resistance 
parameter ruQ [145] . It has been shown for several systems including rod- like 
colloids and linear polymers without self- avoidance [70] that close to equilib- 
rium Gxy ^ '7 and {Gxx ^yy) 

~ 7^, so that Tog is independent of Wi. Our 
results for the orientation resistance are presented in Fig. 17 for various func- 
tionalities / and arm lengths Lf. Data for different Lf collapse onto universal 
curves, which approach a plateau for small shear rates, as expected. For larger 
shear rates, Wi ^ 1, a power-law behavior [142] 



tan(2xG) = 2Gxyl{G, 



XX 



Gyy) = TOc/Wi, 



(94) 



TOG(Wi) - /" Wi'", 



(95) 
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is obtained with respect to the Weissenberg number and the functionahty, 
where a = 0.27 ± 0.02 and fj, = 0.65 ± 0.05. For self-avoiding linear polymers, 
a somewhat smaller exponent /i = 0.54 ± 0.03 was obtained in Refs. [87, 
146], whereas theoretical calculations predict ~ Wi^^"^ in the limit of large 
Weissenberg numbers [147]. 

The data for the average orientation and deformation of a star polymer 
described so far seem to indicate that the properties vary smoothly and 
monotonically from linear polymers to star polymers of high functionality. 
However, this picture changes when the dynamical behavior is considered. 
It is well known by now that linear polymers show a tumbling motion in 
flow, with alternating collapsed and stretched configurations during each cy- 
cle [146,148-150]. For large Weissenberg numbers, this leads to very large 
fluctuations of the largest intramolecular distance of a linear polymer with 
time, as demonstrated experimentally in Refs. [146, 150], and reproduced in 
the MPC simulations [142], sec Fig. 18. A similar behavior is found for / = 3. 
However, for / > 5, a quantitatively different behavior is observed as displayed 
in Fig. 19. Now, the fluctuations of the largest intramolecular distance are 
much smaller and decrease with increasing Weissenberg number as shown in 
Fig. 18, and the dynamics resembles much more the continuous tank-treading 
motion of fluid droplets and capsules. The shape and orientation of such stars 
depends very little on time, while the whole object is rotating. On the other 
hand, a single, selected arm resembles qualitatively the behavior of a linear 
polymer — it also collapses and stretches during the tank-treading motion. The 
successive snapshots of Fig. 20 illustrate the tank-treating motion. Following 
the top left red polymer in the top left image, we see that the extended poly- 
mer collapses in the course of time, moves to the right and stretches again. In 
parallel, other polymers exhibit a similar behavior on the bottom side. More- 
over, the images show that the orientation of a star hardly changes in the 
course of time. 

The rotational dynamics of a star polymer can be characterized quantita- 
tively by calculating the rotation frequency 

z 

^o. = Y.^O-lL0) (96) 

of a star, where 

i=l 

is the instantaneous moment-of-inertia tensor and Ljj is the instantaneous an- 
gular momentum. Since the rotation frequency for all kinds of soft objects — 
such as rods, linear polymers, droplets and capsules — depends linearly on 7 
for small shear rates, the reduced rotation frequency a;/7 is shown in Fig. 21 
as a function of the Weissenberg number. The data approach a;/7 = 1/2 
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Fig. 18. Temporal evolution of the largest intramolecular distance e = max^j |ri — rj] 
of a linear polymer and a star polymer with 50 arms, for the Weissenberg numbers 
Wi = 0.64 and Wi — 64. In both cases, L/ = 40. The time t is measured in units of 
the colhsion time St. corresponds to the fully stretched arms. From Ref. [142]. 
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Fig. 19. Widths of the distribution functions of the largest intramolecular distances, 
ae ~ {{e^) — {e)^)/{e)^, of a Hnear polymer and star polymers with up to 50 arms 
as a function of the Weissenberg number. From Ref. [142] . 

for small Wi, as expected [87, 151]. For larger shear rates, the reduced fre- 
quency decreases due to the deformation and alignment of the polymers in 
the flow field. With increasing arm number, the decrease of w/'y at a given 
Weissenberg number becomes smaller, since the deviation from the spherical 
shape decreases. Remarkably, the frequency curves for all stars with / > 5 
are found to collapse onto a universal scaling function when 07/7 is plotted 
as a function of a rescaled Weissenberg number, see Fig. 21. For high shear 
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Fig. 21. Scaled rotation frequencies as function of a rescaled Weissenberg number 
for various functionalities. Dashed and full lines correspond to w/7 = 1/2 for small 
Wi, and ~ 1/Wi for large Wi, respectively. The inset shows the dependence of 
the rescaling factor ip on the functionality. From Ref. [142]. 



rates, 1.^/7 decays as Wi^ , which implies that the rotation frequency becomes 
independent of 7. A similar behavior has been observed for capsules at high 
shear rates [152]. 

The presented results show that star polymers in shear flow show a very 
rich structural and dynamical behavior. With increasing functionality, stars 
in flow change from linear-polymer-like to capsule- like behavior. These macro- 
molecules are therefore interesting candidates to tune the viscoelastic proper- 
ties of complex fluids. 
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11 

Vesicles and Cells in Hydrodynamic Flows 
11.1 

Introduction 

The flow behavior of fluid droplets, capsules, vesicles and cells is of enormous 
importance in science and technology. For example, the coalescence and break- 
up of fluid droplets is essential for emulsion formation and stability. Capsules 
and vesicles are discussed and used as drug carriers. Red blood cells (RBC) 
flow in the blood stream, and may c;oagulate or be torn apart under unfa- 
vorable flow conditions. Red blood cells also have to squeeze through narrow 
capillaries to deliver their oxygen cargo. White blood cells in capillary flow 
adhere to, roll along and detach again from the walls of blood vessels under 
normal physiological conditions; in inflamed tissue, leukocyte rolling leads to 
flrm adhesion and induces an immunological response. 

These and many other applications have induced an intensive theoretical 
and simulation activity to understand and predict the behavior of such soft, 
deformable objects in flow. In fact, there are some general, quahtative proper- 
ties in simple shear flow, which are shared by droplets, capsules, vesicles and 
cells. When the internal viscosity is low and when they are highly deformable, 
a tank-treading motion is observed (in the case of droplets for not too high 
shear rates), where the shape and orientation are stationary, but particles lo- 
calized at the interface or attached to the membrane orbit around the center 
of mass with a rotation axis in the vorticity direction. On the other hand, 
for high internal viscosity or small deformability, the whole object performs a 
tumbling motion, very much like a colloidal rod in shear flow. However, if we 
take a more careful look, then the behavior of droplets, capsules, vesicles and 
cells is quite different. For example, droplets can break up easily at higher 
shear rates, because their shape is determined by the interfacial tension; fluid 
vesicles can deform much more easily then capsules, since their membrane has 
no shear elasticity; etc. We focus here on the behavior of fluid vesicles and red 
blood cells. 

11.2 

Modeling Membranes 
11.2.1 

Modeling Lipid-Bilayer Membranes 

The modeling of lipid bilayer membranes depends very much on the length 
scale of interest. The structure of the bilayer itself or the embedding of mem- 
brane proteins in a bilayer are best studicxl with atomistic models of both 
lipid and water molecules. Molecular dynamics simulations of such models 
are restricted to about 10^ lipid molecules. For larger system sizes, coarse- 
grained models are required [153-155]. Here, the hydrocarbon chains of lipid 
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molecules are described by short polymer chains of Lennard-Jones particles, 
which have a repulsive interaction with the lipid head groups as well as with 
the water molecules, which arc also modeled as single Lcnnard- Jones spheres. 
Very similar models, with Lennard-Jones interactions replaced by linear "soft" 
dissipative-particle dynamics (DPD) potentials, have also been employed in- 
tensively [156 160]. For the investigation of shapes and thermal fluctuations 
of single- or multi-component membranes, the hydrodynamics of the solvent is 
irrelevant. In this case, it can be advantageous to use a solvent-free membrane 
model, in which the hydrophobic effect of the water molecules is replaced by an 
effective attraction among the hydrocarbon chains [161-164]. This approach is 
advantageous in the case of membranes in dilute solution, because it reduces 
the number of molecules and thus the degrees of freedom to be simulated 
by orders of magnitude. However, it should be noticed that the basic length 
scale of atomistic and coarse-grained or solvent-free models is still on the same 
order of magnitude. 

In order to simulate larger systems, such as giant unilamellar vesicles 
(GUV) or red blood cells (RBC), which have a radius on the order of sev- 
eral micrometers, a different approach is required. It has been shown that in 
this limit the properties of lipid bilayer membranes are described very well by 
modeling the membrane as a two-dimensional manifold embedded in three- 
dimensional space, with the shape and fluctuations controlled by the curvature 
elasticity [165], 



where H = {ci + C2) /2 is the mean curvature, with the local principal curva- 
tures Ci and C2, and the integral is over the whole membrane area. To make 
the curvature elasticity amenable to computer simulations, it has to be dis- 
cretized. This can be done either by using triangulated surfaces [166, 167], or 
by employing particles with properly designed interactions which favor the 
formation of self-assembled, nearly planar sheets [168,169]. In the latter case, 
both scalar particles with isotropic multi-particle interactions (and a curvature 
energy obtained from a moving least-squares method) [169] as well as particles 
with an internal spin variable and anisotropic, multi-body forces [168] have 
been employed and investigated. 

11.2.2 

Dynamically Triangulated Surfaces 

In a dynamically-triangulated surface model [166,167,170-172] of vesicles and 
cells, the membrane is described by Nmb vertices which are connected by 
tethers to form a triangular network of spherical topology, see Fig. 22. The 
vertices have excluded volume and mass nimb- Two vertices connected by 
a bond have an attractive interaction, which keeps their distance below a 
maximum separation Iq. A short-range repulsive interaction among all vertices 




(98) 



Multi-Particle Collision Dynamics 



65 



makes the network self-avoiding and prevents very short bond lengths. The 
curvature energy can be discretized in different ways [f66,f73]. In particular, 
the discretization [173, 174] 



Uc^^yAiE"—"^ (99) 



has been found to give reliable results in comparison with the continuum 
expression (98). Here, the sum over is over the neighbors of a vertex 
i which are connected by tethers. The bond vector between the vertices i 
and j is r^j, and ri,j = \vij\. The length of a bond in the dual lattice is 
Cj.j = ''i,i[cot(6'i) + cot(6'2)]/2, where the angles 9i and 02 are opposite to 
bond ij in the two triangles sharing this bond. Finally, cr^ = 0.25 ^^^^-j (Tijri j 
is the area of the dual cell of vertex i. 



Fig. 22. Triangulated-network model of a fluctuating membrane. AH vertices have 
short-range repulsive interactions symbolized by hard spheres. Bonds represent at- 
tractive interactions with imply a maximum separation £q of connected vertices. 
From Ref. [175]. 



To model the fluidity of the membrane, tethers can be flipped between 
the two possible diagonals of two adjacent triangles. A number ipNt of bond- 
flip attempts is performed with the Metropolis Monte Carlo method [173] at 
time intervals AtsF, where = 3{Nmb — 2) is the number of bonds in the 
network, and < i/^ < 1 is a parameter of the model. Simulation results show 
that the vertices of a dynamically triangulated membrane show diffusion, i.e., 
the squared distance of two initially neighboring vertices increases linearly in 
time. 
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11.2.3 

Vesicle Shapes 

Since the solubility of lipids in water is very low, the number of lipid molecules 
ill a membrane is essentially constant over typical experimental time scales. 
Also, the osmotic pressure generated by a small number of ions or macro- 
molecules in solution, which cannot penetrate the lipid bilaycr, keeps the in- 
ternal volume essentially constant. The shape of fluid vesicles [176] is therefore 
determined by the competition of the curvature elasticity of the membrane, 
and the constraints of constant volume V and constant surface area S. In the 
simplest case of vanishing spontaneous cnirvature. the cnirvature elasticity is 
given by Eq. (98). In this case, the vesicle shape in the absence of thermal 
fluctuations depends on a single dimensionless parameter, the reduced volume 
V* = V/Vo, where Vq = (47r/3)i?g and Ro = (S'/47r)i/2 are the volume and 
radius of a sphere of the same surface area S, respectively. The calculated vesi- 
cle shapes are shown in Fig. 23. There are three phases. For reduced volumes 
not too far from the sphere, elongated prolate shapes are stable. In a small 
range of reduced volumes of V* e [0.592, 0.651], oblate discocyte shapes have 
the lowest curvature energy. Finally, at very low reduced volumes, cup-like 
stomatocyte shapes axe found. 

These shapes are very well reproduced in simulations with dynamically 
triangulated surfaces [177-180]. 



Fig. 23. Shapes of fluid vesicles as a function of the reduced volume V*. D and 

Z)"*" denote the discontinuous prolate-oblate and oblate-stomatocyte transitions, 
respectively. All shapes display rotational symmetry with respect to the vertical 
axis. Prom Ref. [176]. 



11.2.4 

Modeling Red Blood Cells 

Red blood cells have a biconcave disc shape, which can hardly be distinguished 
from the discocyte shape of fluid vesicles with reduced volume V* ~ 0.6, 
compare Fig. 23. However, the membrane of red blood cells is more complex. 
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since a spectrin network is attached to the plasma membrane [181], which 
helps to retain the integrity of the cell in strong shear gradients or capillary 
flow. Due to the spectrin network, the red blood cell membrane has a non-zero 
shear modulus yit. 

The bending rigidity k of RBCs has been measured by micropipette aspira- 
tion [182] and atomic force microscopy [183] to be approximately k = bOksT. 
The shear modulus of the composite membrane, which is induced by the 
spectrin network, has been determined by several techniques; it is found to be 
H = 2 X 10^^ N/m from optical tweezers manipulation [184], while the value 
/z = 6 X 10^^ N/m is obtained from micropipette aspiration [182]. Thus, the 
dimensionless ratio IjlRq/k ~ 100, which implies that bending and stretching 
energies arc roughly of equal importance. 

Theoretically, the shapes of RBCs in the absence of flow have been cal- 
culated very successfully on the basis of a mechanical model of membranes, 
which includes both curvature and shear elasticity [185,186]. In particular, 
it has been shown recently that the full stomatocyte-discocyte-echinocyte se- 
quence of RBCs can be reproduced by this model [186]. 

The composite membrane of a red blood cell, consisting of the lipid bilayer 
and the spectrin network, can be modeled as a composite network, which 
consists of a dynamically-triangulated surface as in the case of fluid vesicles, 
coupled to an additional network of harmonic springs with fixed connectivity 
(no bond-flip) [185, 187]. Ideally, the bond length of the elastic network is 
larger than of the fluid mesh [185] — in order to mimic the situation of the red 
blood cell membrane, where the average distance of anchoring points is about 
70 nm, much larger than, the size of a lipid molecule — and thereby allow, 
for example, for thermal fluctuations of the distances between neighboring 
anchoring points. On the other hand, to investigate the behavior of cells on 
length scales much larger the mesh size of the spectrin network, it is more 
eflicient to use the same number of the bonds for both the fluid and the 
tethered networks [187]. The simplest case is a harmonic tethering potential, 
(l/2)/eei(ri — Vj)"^. This tether network generates a shear modulus ji = v^fcei- 

11.3 

Modeling Membrane Hydrodynamics 

Solvent-free models, triangulated surfaces and other discrctizcd curvature 
models have the disadvantage that they do not contain a solvent, and therefore 
do not describe the hydrodynamic behavior correctly. However, this apparent 
disadvantage can be turned into an advantage by combining these models with 
a mesoscopic hydrodynamics technique. This approach has been employed for 
dynamically triangulated surfaces [37, 180] and for meshless membrane models 
in combination with MPC [188], as well as for fixed membrane triangulations 
in combination with both MPC [187] and the Lattice Boltzmann method 
(LBM) [189]. 
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The solvent particles of the MPC fluid interact with the membrane in two 
ways to obtain an impermeable membrane with no-slip boundary conditions. 
First, the membrane vertices arc included in the MPC collision procedure, 
as suggested for polymers in Ref. [68], compare Sec. 7.1. Second, the solvent 
particles are scattered with a bounce-back rule from the membrane surface. 
Here, solvent particles inside (1 < 2 < A^m) and outside {Nin < i < Ng) 
of the vesicle have to be distinguished. The membrane triangles are assumed 
to have a finite but very small thickness 6 = 2lhs- The scattering process 
is then performed at discrete time steps AIbs, so that scattering docs not 
occur exactly on the membrane surface, but the solvent particles can penetrate 
slightly into the membrane film [180] . A similar procedure has bee suggested in 
Ref. [26] for spherical colloidal particles embedded in a MPC solvent. Particles 
which enter the membrane film, i.e., which have a distance to the triangulated 
surface smaller than lf,s, or interior particles which reach the exterior volume 
and vice versa, are scattered at the membrane triangle with the closest center 
of mass. Explicitly [180], 

V^-)(t) = Mt) - — ^^^(V.W - Vtriit)) (100) 

ms + 6mmb 

when {\s{t) — 'Vtri{t)) ■ ntri < 0, where Vs(f) and "Vtriit) are the velocities 
of the solvent particle and of the center of mass of the membrane triangle, 
respectively, and nj,.j is the normal vector of the triangle, which is oriented 
towards the outside (inside) for external (internal) particles. 

The bond flips provide a very convenient way to vary the membrane vis- 
cosity rirnb, which increases with decreasing probability ijj for the selection 
of a bond for a bond-flip attempt. The membrane viscosity has been deter- 
mined quantitatively from a simulation of a flat membrane in two-dimensional 
Poiseuille flow. The triangulated membrane is put in a rectangular box of size 

X Ly with periodic boundary conditions in the x-direction. The edge ver- 
tices at the lower and upper boundary {y = ±Ly/2) are flxed at their positions. 
A gravitational force {m,rnb9, 0) is applied to all membrane vertices to induce 
a flow. Rescaling of relative velocities is employed as a thermostat. Then, 
the membrane viscosity is calculated from r]mb = PmbdLy/Svmax, where pmb 
is average mass density of the membrane particles, and v„iax the maximum 
velocity of the parabolic flow proflle. The membrane viscosity r/mfa, which is 
obtained in this way, is shown in Fig. 24. As tjj decreases, it takes longer and 
longer for a membrane particle to escape from the cage of its neighbors, and 
rjmb increases. This is very similar to the behavior of a hard-sphere fluid with 
increasing density. Finally, for = 0, the membrane becomes solid. 
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Fig. 24. Dependence of the membrane viscosity ijmb on the probability ip for the 
selection of a bond for a bond-flip attempt, for a membrane with Nmb = I860 
vertices. EYom Ref. [180]. 

11.4 

Fluid Vesicles in Shear Flow 

The dynamical behavior of fluid vesicles in simple shear flow has been stud- 
ied experimentally [190-193], theoretically [194-201], numerically with the 
boundary-integral technique [202,203] or the phase-field method [203,204], 
and with mesoscale solvents [37,180,205]. The vesicle shape is now deter- 
mined by the competition of the curvature elasticity of the membrane, the 
constraints of constant volume V and constant surface area S, and the exter- 
nal hydrodynamic forces. 

Shear flow is characterized (in the absence of vesicles or cells) by the flow 
field V = jyex, where e^; is a unit vector, compare Sec. 10.4. The control 
parameter of shear flow is the shear rate 7, which has the dimension of an 
inverse time. Thus, a dimensionless, scaled shear rate 7* = -yr can be deflned, 
where r is the longest relaxation time of a vesicle. It is given by r = tjoRq/k, 
where rjQ is the solvent viscosity, Rq the average radius, and k the bending 
rigidity [206]. For 7* < 1, the internal vesicle dynamics is fast compared to 
the external perturbation, so that the vesicle shape is hardly affected by the 
flow field, whereas for 7* > 1, the flow forces dominate and the vesicle is in a 
non-equilibrium steady state. 

One of the difficulties in theoretical studies of the hydrodynamic effects 
on vesicle dynamics is the no-slip boundary condition for the embedding fluid 
on the vesicle surface, which changes its shape dynamically under the effect 
of flow and curvature forces. In early studies, a fluid vesicle was therefore 
modeled as an ellipsoid with fixed shape [194]. This simplified model is still 
very useful as a reference for the interpretation of simulation results. 
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11.4.1 

Generalized Keller-Skalak Theory 

The theory of Keller and Skalak [194] describes the hydrodynamic behavior of 
vesicles of fixed ellipsoidal shape in shear flow, with the viscosities rjin and rjo of 
the internal and external fluids, respectively. Despite of the approximations 
needed to derive the equation of motion for the inclination angle 6, which 
measures the deviation of the symmetry axis of the ellipsoid with the flow 
direction, this theory describes vesicles in flow surprisingly well. It has been 
generalized later [197] to describe the effects of a membrane viscosity rjmb- 

The main result of the theory of Keller and Skalak is the equation of 
motion for the inclination angle [194], 

|^=i^{-l + Bcos(2^)}, (102) 

where i? is a constant, which depends on the geometrical parameters of the 
ellipsoid, on the viscosity contrast r]*^ = r]in/r]o, and the scaled membrane 
viscosity rj^^ = r]mb/{iloRo) [180,194,197], 

where /o, /i, /2, and are geometry- dependent parameters. In the spherical 
limit, B ^ CO. Eq. (102) implies the following behavior: 

• For B > \, there is a stationary solution, with cos(20) = 1/B. This cor- 
responds to a tank-treading motion, in which the orientation of the vesicle 
axis is time independent, but the membrane itself rotates around the vor- 
ticity axis. 

• For < 1, no stationary solution exists, and the vesicle shows a tumbling 
motion, very similar to a solid rod-like colloidal particle in shear flow. 

• The shear rate 7 only determines the time scale, but does not affect the 
tank-treading or tumbling behavior. Therefore, a transition between these 
two types of motion can only be induced by a variation of the vesicle shape 
or the viscosities. 

However, the vesicle shape in shear flow is often not as constant as as- 
sumed by Keller and Skalak. In these situations, it is very helpful to compare 
simulation results with a generalized Keller-Skalak theory, in which shape de- 
formation and thermal fluctuations are taken into account. Therefore, a phe- 
nomenological model has been suggested in Ref. [180] , in which in addition to 
the inclination angle 6 a sec;ond parameter is introduced to characterize the 
vesicle shape and deformation, the asphericity [207] 

^ ^ (Ai-A2)^ + (A2-A3)^ + (A3-Ai)^ ^ ^^^^^ 

where Ai < A2 < A3 are the eigenvalues of the moment-of-inertia tensor and 
the squared radius of gyration iJ^ = Ai -|- A2 -|- A3. This implies a = for 
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spheres (with Ai = A2 = A3), a = 1 for long rods (with Ai = A2 <C A3), and 
a = 1/4 for flat disks (with Ai ^ A2 = A3). The generalized Keller-Skalak 
model is then defined by the stochastic equations 

Ca^a = -dF/da + A7sin(2^) + CaQait) (105) 

±e = ^7{-l + B{a) cos{2e)} + ge{t), (106) 
with Gaussian white noises ga and gg, which are determined by 

{9a{t)) - {ge{t)) = (.ga(t),%(t')> = 0, 

{ga{t)gc{t')) = 2D^S{t-t'), (107) 
{g9{t)ge{t'))^2De6{t-t'), 

friction coefficients (a and (g, and diffusion constants = knT/^a and 
Dg = k^T/Q). Note that does not appear in Eq. (106); it drops out because 
the shear force is also caused by friction. 

The form of the stochastic equations (105) and (106) is motivated by the 
following considerations. The first term in Eq. (105), dF/da, is the thermo- 
dynamic force due to bending energy and volume constraints; it is calculated 
from the free energy F{a). The second term of Eq. (105) is the deformation 
force due to the shear flow. Since the hych'odynamic forces elongate the vesi- 
cle for < ^ < tt/2 but push to reduce the elongation for — 7r/2 < 9 < 0, 
the flow forces should be proportional to sin(2^) to leading order. The ampli- 
tude A is assumed to be independent of the asphericity a. (a and A can be 
estimated [205] from the results of a perturbation theory [199] in the quasi- 
spherical limit. Eq. (106) is adapted from Keller-Skalak theory. While B is a 
constant in Keller-Skalak theory, it is now a function of the (time-dependent) 
asphericity a in Eq. (106). 



11.4.2 

Effects of Membrane Viscosity: Tank- Treading and Tumbling 

The theory of Keller and Skalak [194] predicts for fluid vesicles a transition 
from tank-treading to tumbling with increasing viscosity contrast rjin/rja. This 
has been confirmed in recent simulations based on a phase- field model [203]. 

The membrane viscosity rjmb is also an important factor for the vesicle dy- 
namics in shear fiow. For example, the membrane of red blood cells becomes 
more viscous on aging [197,208] or in diabetes mellitus [209]. Experiments 
indicate that the energy dissipation in the membrane is larger than that in- 
side a red blood cell [196, 197]. Furthermore, it has been shown recently that 
vesicles can not only be made from lipid bilayers, but also from bilayers of 
block copolymers [210]. The membrane viscosity of these "polymersomes" is 
several orders of magnitude larger than for liposomes, and can be changed 
over a wide range by varying the polymer chain length [211]. 
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Fig. 25. Snapshot of a discocyte vesicle in shear flow with reduced shear rate 
7* = 0.92, reduced volume V* = 0.59, membrane viscosity rj^i, — and viscosity 
contrast rji„/rjo — 1. The arrows represent the velocity field in the xz-plane. From 
Ref. [180]. 

A variation of the membrane viscosity can be implemented easily in dy- 
namically triangulated surface models of membranes, as explained in Sec. 11.2.2 
An example of a discocyte in tank-treading motion, which is obtained by such 
a membrane model [180], is shown in Fig. 25. Simulation results for the inclina- 
tion angle as a function of the reduced membrane viscosity Ty^j, = 77m6/('7o-Ro) 
are shown in Fig. 26. This demonstrates the tank-treading to tumbling transi- 
tion of fluid vesicles with increasing membrane viscosity. The threshold shear 
rate decreases with decreasing reduced volume V* , since with increasing de- 
viation from the spherical shape, the energy dissipation within the membrane 
increases. Interestingly, the discocyte shape is less affected by the membrane 
viscosity than the prolate shape for V* — 0.59, since it is more compact — in 
contrast to a vesicle with viscosity contrast rjm/vo > 1- where the prolate 
shape is affected less [180]. 

Figure 26 also shows a comparison of the simulation data with results of 
K-S theory for fixed shape, both without and with thermal fluctuations. Note 
that there are no adjustable parameters. The agreement of the results of the- 
ory and simulations is excellent in the case of vanishing membrane viscosity, 
Vmb = 0. For small reduced volumes, V* ~ 0.6, the tank-treading to tum- 
bling transition is smeared out by thermal fluctuations, with an intermittent 
tumbling motion occurring in the crossover region. This behavior is captured 
very well by the generalized K-S model with thermal fluctuations. For larger 
reduced volumes and non-zero membrane viscosity, signiflcant deviations of 
theory and simulations become visible. The inclination angle 9 is found to de- 
crease much more slowly with increasing membrane viscosity than expected 
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Fig. 26. Average inclination angle (9) as a function of reduced membrane viscosity 
ri^tj for the shear rate -y* — 0.92 and various reduced volumes V* . Results are 
presented for prolate (circles) and discocyte (squares) vesicles with V* = 0.59, as 
well as prolate vesicles with V* = 0.66 (triangles), 0.78 (diamonds), 0.91 (crosses), 
and V — 0.96 (pluses). The solid and dashed lines are calculated by K-S theory, 
Eqs. (102) and (103), for prolate {V* = 0.59, 0.66, 0.78, 0.91, and 0.96) and oblate 
ellipsoids {V* — 0.59), respectively. The dashed-dotted lines are calculated from 
Eq. (106) with thermal fluctuations, for V* = 0.66, V* = 0.78, and the rotational 
Peclet number "//Dg — 600 (where Dg is the rotational diffusion constant). From 
Ref. [180]. 

theoretically. This is most likely due to thermal membrane undulations, which 
are not taken into account in K-S theory. 
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Fig. 27. Temporal evolution of vesicle deformation an and inclination angle 0, for 
V =0.78 and 77*jf, = 2.9. Here, cto = {L\— L2)/ {Li +L2), where Li and L2 are the 
maximum lengths in the direction of the eigenvectors of the gyration tensor in the 
vorticity plane. The solid (red) and dashed (blue) lines represent simulation data 
for 7* = 3.68 and 0.92 (K/fcsT = 10 and 40, with 7770^0/^13^ = 36.8), respectively. 
The solid lines in (b) are obtained from Eqs. (105), (106) without thermal noise for 
7* = 1.8, 3.0, and 10 (from top to bottom). From Ref. [205]. 



11.4.3 

Swinging of Fluid Vesicles 

Recently, of a new type of vesicle dynamics in shear flow has been observed 
experimentally [192], which is characterized by oscillations of the inclination 
angle 9 with 9{t) E [—9o,6o] and 9o < tt/2. The vesicles were found to 
transit from tumbling to this oscillatory motion with increasing shear rate 
7. Simultaneously with the experiment, a "vacillating- breathing" mode for 
quasi-spherical fluid vesicles was predicted theoretically, based on a spherical- 
harmonics expansion of the equations of motion to leading order (without 
thermal fluctuations) [199]. This mode exhibits similar dynamical behavior 
as seen experimentally; however, it "coexists" with the tumbling mode, and 
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Fig. 28. Dynamical phase diagram as a function of viscosity contrast ?7*„ — rjin/rio, 
for ry^ij = and various reduced volumes V* , calculated from Eqs. (105), (106) 
without thermal noise. The tank-treading phase is located on the left-hand-side 
of the dashed lines. The solid lines represent the tumbling-to-swinging transitions. 
From Ref. [205]. 



its orbit depends on the initial deformation, i.e., it is not a limit cycle. Fur- 
thermore, the shear rate appears only as the basic time scale, and therefore 
cannot induce any shape transitions. Hence it does not explain the tumbling- 
to-oscillatory transition seen in the experiments [192]. 

Simulation data for the oscillatory mode — which has also been denoted 
"trembling" [192] or "swinging" [205] mode — are shown in Fig. 27. The sim- 
ulation results demonstrate that the transition can indeed be induced by in- 
creasing shear rate, and that it is robust to thermal fluctuations. Figure 27 
also shows that the simulation data are well captured by the generalized K- 
S model, Eqs. (105) and (106), which takes into account higher-order con- 
tributions in the curvature energy of a vesicle. The theoretical model can 
therefore be used to predict the full dynamic phase diagram of prolate vesi- 
cles as a function of shear rate and membrane viscosity or viscosity contrast, 
compare Fig. 28. The swinging phase appears at the boundary between the 
tank-treading and the tumbling phase for sufficiently large shear rates. The 
phase diagram explains under which conditions the swinging phase can be 
reached from the tumbling phase with increasing shear rate — as observed ex- 
perimentally [192]. 

The generalized K-S model is designed to capture the vesicle flow behavior 
for non-spherical shapes sufficiently far from a sphere. For quasi-spherical 
vesicles, a derivation of the equations of motion by a systematic expansion 
in the undulation amplitudes gives quantitatively more reliable results. An 
expansion to third order results in a phase diagram [200,201], which agrees 
very well with Fig. 28. 



76 



G. Gompper, T. Ihle, D.M. KroU, and R.G. Winkler 



T — I — I — I — I — I — r~?i — I — I — I — I — I — I — I — r 
+ >■ X > 

transformation from 
/ Etomatocyte to prolate _ 

^ ^- 

.*L h ' tanktreading-tumbling 

/ transition o_ 

. — - — transformation from- 

I discocyte to prolate 

1(p~ I ° ° 



J I I U I I I I I I I I I L 



0,5 * 1 1.5 



Fig. 29. Dynamical phase diagram of a vesicle in shear flow for reduced volume 
V — 0.59. Symbols show simulated parameter values, and indicate tank-treading 
discocyte and tank-treading prolate (o), tank-treading prolate and unstable disco- 
cyte (A), tank-treading discocyte and tumbling (transient) prolate (□), tumbling 
with shape oscillation (o), unstable stomatocyte stable stomatocyte (x), and 
near transition (■). The dashed lines are guides to the eye. From Ref. [180]. 



11.4.4 

Flow-Induced Shape Transformations 

Sliear flow does not only induce different dynamical modes of prolate and 
oblate fluid vesicles, it can also induce phase transformations. The simplest 
case is a oblate fluid vesicle with rj^^b — and viscosity contrast 77m/f?o = 1- 
When the reduced shear rate reaches 7* ~ 1, the discocyte vesicles are 
stretched by the flow forces into a prolate shape [37, 180,202]. A similar tran- 
sition is found for stomatocyte vesicles, except that in this case a larger shear 
rate 7* ~ 3 is required. In the case of non-zero membrane viscosity, a rich 
phase behavior appears, see Fig. 29. 

Surprisingly, flow forces can not only stretch vesicles into a more elongated 
shape, but can also induce a transition from an elongated prolate shape into 
a more compact discocyte shape [180]. Simulation results for the latter tran- 
sition are shown in Fig. 30. This transition is possible, because in a range of 
membrane viscosities, the prolate shape is in the tumbling phase, while the 
oblate shape is tank-treading, compare Fig. 26. Of course, this requires that 
the free energies of the two shapes are nearly equal, which implies a reduced 
volume of V* ~ 0.6. Thus, a prolate vesicle in this regime starts tumbling; 
as the inclination angle becomes negative, shear forces push to shrink the 
long axis of the vesicle; when this force is strong enough to overcome the 
free-energy barrier between the prolate and the oblate phase, a shape trans- 
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formation can be induced, compare Fig. 30. The vesicle then remains in the 
stable tank-treading state. 

At higher shear rates, on intermittent behavior has been observed, in which 
the vesicle motion changes in irregular intervals between tumbling and tank- 
treading [180]. 
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Fig. 30. Time dependence of asphericity a and inclination angle 6, for 7* — 1.84, 
'7mb ~ 1-62, and V* — 0.59. The dashed lines are obtained from Eqs. (105) and 
(106), with Co = 100, A = 12, and B{a) = 1.1 - 0.17a. From Ref. [212]. 



11.4.5 

Vesicle Fluctuations under Flow in Two Dimensions 

At finite temperature, stochastic fluctuations of the membrane due to ther- 
mal motion affect the dynamics of vesicles. Since the calculation of thermal 
fluctuations under flow conditions requires long times and large membrane 
sizes (in order to have a sufficient range of undulation wave vectors), simu- 
lations have been performed for a two-dimensional system in the stationary 
tank-treading state [213]. For comparison, in the limit of small deformations 
from a circle, Langevin-type equations of motion have been derived, which are 
highly nonlinear due to the constraint of constant perimeter length [213]. 

The effect of the shear flow is to induce a tension in the membrane, which 
reduces the amplitude of thermal membrane undulations. This tension can 
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be extracted directly from simulation data for the undulation spectrum. The 
reduction of the undulation amplitudes also implies that the fluctuations of 
the inclination angle 9 get rc(hic;c(l with increasing shear rate. The theory 
for quasi-circular shapes predicts a universal behavior as a function of the 
scaled shear rate 7*Zi^/^K/ {RokBT), where 7* = 7?7o-Ro/'^ reduced shear 
rate in two dimensions, and A is the dimcnsionlcss excess membrane length. 
Theory and simulation results for the inclination angle as a function of the 
reduced shear rate are shown in Fig. 31. There are no adjustable parameters. 
The agreement is excellent as long as the deviations from the circular shape 
are not too large [213]. 



A 

N 

CD 
< 
V 




0.1- 



7 A K/R^kj^T 



Fig. 31. Fluctuations of the inclination angle (AO^y^'^ of a two-dimensional fluc- 
tuating vesicle in shear flow, as a function of scaled shear rate 7*Zi^^'^K/(_Rofcsr), 
where A is the dimensionless excess membrane length. Symbols indicate simula- 
tion data for different internal vesicle areas A for fixed membrane length, with 
A* = A/nRo = 0.95 (squares), A* = 0.90 (triangles). A* = 0.85 (stars), and 
A* = 0.7 (circles). The solid line is the theoretical result in the quasi-circular limit. 
Prom Ref. [213]. 
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11.5 



Fluid Vesicles and Red Blood Cells In Capillary Flow 
11.5.1 

RBC Deformation In Narrow Capillaries 

The deformation of single RBCs and single fluid vesicles in capillary flows were 

studied theoretically by lubrication theories [214-216] and boundary-integral 
methods [217-219]. In most of these studies, axisymmetric shapes which are 
coaxial with the center of the capillary were assumed and cylindrical coordi- 
nates were employed. In order to investigate non- axisymmetric shapes as well 
as flow-induced shape transformations, a fully three-dimensional simulation 
approach is required. 

We focus here on the behavior of single red blood cells in capillary flow 
[187], as described by a triangulated surface model for the membrane (compare 
Sec. 11.2.4) immersed in a MFC solvent (see Sec. 11.2.4). The radius of the 
capillary, Reap, is taken to be slightly larger than the mean vesicle or RBC 
radius, Rq = \/Sj4^, where S is the membrane area. Snapshots of vesicle 
and RBC shapes in flow are shown in Fig. 32 for a reduced volume of V* = 
0.59, where the vesicle shape at rest is a discocytc. For sufficiently small flow 
velocities, the discocyte shape is retained. However, the discocyte is found not 
in a coaxial orientation; instead the shortest eigenvalue of the gyration tensor 
is oriented perpendicular to the cylinder axis [187]. Since two opposite sides of 
the rim of the discocyte are closer to wall where the flow velocity is small, the 
rotational symmetry is slightly disturbed and the top view looks somewhat 
triangular, sec Fig. 32A. With increasing flow velocity, a shape transition to an 
axisymmetric shape occurs. In the case of fluid vesicles this is a prolate shape, 
while in the case of RBCs a parachute shape is found. Such parachute shapes 
of red blood cells have previously been observed experimentally [209,220]. 

The fundamental difference between the flow behaviors of fluid vesicles 
and red blood cells at high flow velocities is due to the shear elasticity of the 
spectrin network. Its main effect for jiR^/ k > 1 is to suppress the discocytc-to- 
prolate transition, because the prolate shape would acquire an elastic energy 
of order jiR^. In comparison, the shear stress in the parachute shape is much 
smaller. 

Some diseases, such as diabetes mellitus and sickle cell anemia, change the 
mechanical properties of RBCs; a reduction of RBC deformabihty was found 
to be responsible for an enhanced flow resistance of blood [221]. Therefore, it 
is very important to understand the relation of RBC elasticity and flow prop- 
erties in capillaries. The flow velocity at the discocyte-to-prolate transition of 
fluid vesicles and at the discocyte-to-parachute transition is therefore shown 
in Fig. 33 as a function of the bending rigidity and the shear modulus, respec- 
tively. In both cases, an approximately linear dependence is obtained [187], 



T 



= 0.1 



liRl 



K 



(108) 



R. 



'cap 



80 



G. Gompper, T. Ihle, D.M. KroU, and R.G. Winkler 



■A 




side view 



Fig. 32. Snapshots of vesicles in capillary flow, with bending rigidity K/fesT = 
20 and capillary radius Reap = IARq. (A) Fluid vesicle with discoidal shape at 
the mean fluid velocity VmT/Rcap = 41, both in side and top views. (B) Elastic 
vesicle (RBC model) with parachute shape at VmT / Reap ~ 218 (with shear modulus 
nRl/ksT = 110). The blue arrows represent the velocity field of the solvent. (C) 
Elastic vesicle with slipper-like shape at VmT/Rcap ~ 80 (with ^R^/ksT = 110). 
The inside and outside of the membrane are depicted in red and green, respectively. 
The upper front quarter of the vesicle in (B) and the front half of the vesicle in (C) 
are removed to allow for a look into the interior; the black circles indicate the lines 
where the membrane has been cut in this procedure. Thick black lines indicate the 
walls the cylindrical capillary. From Ref. [187]. 

This result suggests that parachute shapes of RBCs should appear for flow 
velocities larger than = SQQRcap / t — 0.2nim/s under physiological condi- 
tions. This is consistent with the experimental results of Ref. [222], and is in 
the range of micro-circulation in the human body. 

Figure 33 (right) also shows that there is a metastable region, where disco- 
cytes are seen for increasing flow velocity, but parachute shapes for decreasing 
flow velocity. This hysteresis becomes more pronounced with increasing shear 
modulus. It is believed to be due to a suppression of thermal fluctuations with 
increasing y^R^/ksT. 
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Fig. 33. Critical flow velocity Vm of the discocyte-to-parachute transition of elastic 
vesicles and of the discocyte-to-prolate transition of fluid vesicles, as a function of 
the bending rigidity k. From Ref. [187]. 

11.5.2 

Flow in Wider Capillaries 

The flow of many red blood cells in wider capillaries has also been investi- 
gated by several simulation techniques. Discrete fluid-particle simulations — 
an extension of dissipative-particle dynamics (DPD) — in combination with 
bulk-elastic discocyte cells (in contrast to the membrane elasticity of real red 
blood cells) have been employed to investigate the dynamical clustering of red 
blood cells in capillary vessels [223,224]. An immersed flnite-element model — a 
combination of the immersed boundary method for the solvent hydrodynam- 
ics [225] and a finite-element method to describe the membrane elasticity — has 
been developed to study red blood cell aggregation [226]. Finally, it has been 
demonstrated that the Lattice-Boltzmann method for the solvent in combina- 
tion with a triangulated mesh model with curvature and shear elasticity for 
the membrane can be used efficiently to simulate RBC suspensions in wider 
capillaries [189]. 

12 

Viscoelastic Fluids 

One of the unique properties of soft matter is its viscoelastic behavior [13]. 
Due to the long structural relaxation times, the internal degrees of freedom 
cannot relax sufficiently fast in an oscillatory shear flow already at moderate 
frequencies, so that there is an elastic restoring force which pushes the system 
back to its previous state. Well studied examples of viscoelastic fluids are 
polymer solutions and polymer melts [13,70]. 
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The viscoelastic behavior of polymer solutions leads to many unusual flow 
phenomena, such as viscoelastic phase separation [227] . There is also a second 
level of complexity in soft matter systems, in which a colloidal component is 
dispersed in a solvent, which is itself a complex fluid. Examples are spherical or 
rod-like colloids dispersed in polymer solutions. Shear flow can induce particle 
aggregation and alignment in these systems [228]. 

It is therefore desirable to generalize the MPC technique to model vis- 
coelastic fluids, while retaining as much as possible of the computational sim- 
plicity of standard MPC for Newtonian fluids. This can be done by replacing 
the point particles of standard MPC by harmonic dumbbells with spring con- 
stant K [229]. 

As for point particles, the MPC algorithm consists of two steps, streaming 
and collisions. In the streaming step, within a time interval At, the motion of 
all dumbbells is governed by Newton's equations of motion. The center-of-mass 
coordinate of each dumbbell follows a simple ballistic trajectory. The evolution 
of the relative coordinates of dumbbell i, which consists of two monomers at 
positions rji(t) and rj2(t) with velocities Vii{t) and Vi2{t), respectively, is 
determined by the harmonic interaction potential, so that 

rii{t + At) - r,2it + At) = A,{t) cos{ujoAt) + Bi{t) sin(a;o^i) ; (109) 
va{t + At) - Vi2{t + At) = -LOoAi{t) sm{uo At) 

+LJaB,{t) cos{ujQ At) , (110) 

with angular frequency Wq = ^JlKjm. The amplitudes Aj(t) and Bj(t) are 
determined by the initial positions and velocities at time t. The collision step 
is performed for the two point particles constituting a dumbbell in exactly 
the same way as for MPC point-particle fluids. This implies, in particular, 
that the various collision rules of MPC, such as SRD, AT-a or AT-|-a, can 
all be employed also for simulations of viscoelastic solvents, depending on the 
requirements of the system under consideration. Since the streaming step is 
only a little more time consuming and the collision step is identical, simula- 
tions of the viscoelastic MPC fluid can be performed with essentially the same 
efficiency as for the standard point-particle fluid. 

The behavior of harmonic dumbbells in dilute solution has been studied in 
detail analytically [230]. These results can be used to predict the zero-shear 
viscosity ry and the storage and loss moduli, G"(a;) and G"(w) in oscillatory 
shear with frequency w, of the MPC dumbbell fluid. This requires the solvent 
viscosity and diffusion constant of monomers in the solvent. Since the vis- 
coelastic MPC fluid consists of dumbbells only, the natural assumption is to 
employ the viscosity rjMPC and diffusion constant D of an. MPC point-particle 
fluid of the same density. The zero-shear viscosity is then found to be [229] 

V = Vmpc + o ' (111) 

where 
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Fig. 34. (a) Storage G' and (b) loss moduli G", as function of oscillation frequency uo 
on a double-logarithmic scale, for systems of dumbbells with various spring constants 
ranging from K = 0.2 to isT = 1.0. Simulations are performed in two dimensions with 
the SRD collision rule. The wall separation and the collision time are Ly = 10 and 
/At = 0.02, respectively. From Ref. [229]. 
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Similarly, the storage and loss modulus, and the average dumbbell extension, 
are predicted to be [229] 
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(113) 



2 1 + (w/wh) 
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G" = riMPC w + 
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(115) 



Simulation data are shown in Fig. 34, together with the theoretical predica- 
tions (113) and (114). The comparison shows a very good agreement. This 
includes not only the linear and quadratic frequency dependence of G" and 
G' for small w, respectively, but also the leveling off when u) reaches ujh- In 
case of G" , there is quantitative agreement without any adjustable param- 
eters, whereas G' is somewhat overestimated by Eq. (113) for small spring 
constants K. The good agreement of theory and simulations implies that the 
characteristic frequency decreases linearly with decreasing spring constant K 
and mean free path A (since D ~ A). A comparison of other quantities, such 
as the zero-shear viscosity, shows a similar quantitative agreement [229] . 



Conclusions and Outlook 

In the short time since Malevanets and Kapral introduced multi-particle col- 
lision (MPC) dynamics as a particle-based mesoscale simulation technique 
for studying the hydrodynamics of complex fluids, there has been enormous 
progress. It has been shown that kinetic theory can be generalized to cal- 
culate transport coefficients, several collision algorithms have been proposed 
and employed, and the method has been generalized to describe multi-phase 
flows and viscoelastic fluids. The primary applications to date — which include 
studies of the equilibrium dynamics and flow properties of colloids, polymers, 
and vesicles in solution — have dealt with mesoscopic particles embedded in 
a single-component Newtonian solvent. An important advantage of this algo- 
rithm is that it is very straightforward to model the dynamics for the embed- 
ded particles using a hybrid MPC-molecular dynamics simulations approach. 
The results of these studies are in excellent quantitative agreement with both 
theoretical predictions and results obtained using other simulation techniques. 

How will the method develop in the future? This is of course difficult to 
predict. However, it seems clear that there will be two main directions, a fur- 
ther development of the method itself, and its application to new problems in 
Soft Matter hydrodynamics. On the methodological front, there are several 
very recent developments, like angular-momentum conservation, multi-phase 
flows and viscoelastic fluids, which have to be explored in more detail. It 
will also be interesting to combine them to study, e.g., multi-phase flows of 
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viscoelastic fluids. On the application side, the trend will undoubtedly be to- 
wards more complex systems, in which thermal fluctuations are important. In 

such systems, the method can play out its strengths, because the interactions 
of colloids, polymers, and membranes with the mesoscale solvent can all be 
treated on the same basis. 
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